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more conservative point of view to understand the evo- 

The present work aims to push analytical modeling of mtion 01 structure on galaxy cluster scales (for a recent 

structure formation into a regime that may be placed be- theoretical attempt to address halo structure within ki- 

tween the formation epoch of large-scale structure and nctic theory see Ma & Bertschinger 2004). 

the onset of virialization of gravitationally bound ob- The current status of analytical models concern- 

jects. Phenomenologically, this regime is characterized by ing large-scale structure formation may be centred on 

a stabilization of structures that formed out of gravita- Zel'dovich's approximation together with its foundations 

tional instability of a cold medium. The physical origin (the Lagrangian perturbation theory) and optimizations 

of this adhesive clustering effect is the balance of gravita- using filtering techniques for initial perturbation spec- 

tional forces and dynamical stresses in collisionless matter, tra (Zel'dovich 1970, 1973; Shandarin & Zel'dovich 1989; 

The latter arise from a subsequently establishing multi- Buchert 1989, 1992; Melott 1994; Bouchet et al. 1992, 

stream hierarchy within collapsing high-density regions. 1995; Sahni & Coles 1995; Buchert 1996; Ehlers & Buchert 

Although multi-stream forces tend to disperse structures, 1997). These schemes are capable of modeling the evolu- 

the resulting effect together with gravity tends to stabi- tion of generic spectra including Cold-Dark-Matter cos- 

lize them. This regime is summarized by the term non- mogonies down to galaxy cluster scales (Bouchet et al. 

dissipative gravitational turbulence advanced by Gurevich 1995; Buchert et al. 1994, Melott et al. 1995; WeiB et al. 

& Zybin (1995). The models we investigate are a focus of 1997; Hamana 1998). Above these scales the optimized 

current research, since efforts to simulate Hubble volumes Lagrangian schemes roughly reproduce the results of N- 

of the Universe and to understand galaxy halo formation body simulations. Since the exact solution is not known, 
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Abstract. The notion of adhesion has been advanced for the phenomenon of stabilization of large-scale structure 
emerging from gravitational instability of a cold medium. Recently, the physical origin of adhesion has been 
identified: a systematic derivation of the equations of motion for the density and the velocity fields leads naturally 
to the key equation of the 'adhesion approximation' - however, under a set of strongly simplifying assumptions. 
In this work, we provide an evaluation of the current status of adhesive gravitational clustering and a clear 
explanation of the assumptions involved. Furthermore, we propose systematic generalizations with the aim to 
relax some of the simplifying assumptions. We start from the general Newtonian evolution equations for self- 
gravitating particles on an expanding Friedmann background and recover the popular 'dust model' (pressureless 
fluid) , which breaks down after the formation of density singularities; then we investigate, in a unified framework, 
two other models which, under the restrictions referred to above, lead to the 'adhesion approximation'. We 
apply the Eulerian and Lagrangian perturbative expansions to these new models and, finally, we discuss some 
non-perturbative results that may serve as starting points for workable approximations of non-linear structure 
formation in the multi-stream regime. In particular, we propose a new approximation that includes, in limiting 
cases, the standard 'adhesion model' and the Eulerian as well as Lagrangian first-order approximations. 

Key words. Gravitation - Methods: analytical - Cosmology: theory - (Cosmology:) large-scale structure of 
Universe 

1. Gravitational clustering are faced within a single approach. However, we take a 
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we seek agreement between the two modeling techniques 
used. Below these scales, both techniques tend to fail as 
a result of both poorly understood physics and poor res- 
olution power, respectively. Agreement between N-body 
simulations and simple solutions of Lagrangian perturba- 
tion schemes may be considered as supporting analytical 
models, but they could as well be considered as a draw- 
back of N-body simulations given the simplicity of the an- 
alytical approximations compared with the complexity of 
non-linear self-gravity. A deeper understanding of struc- 
ture formation below galaxy cluster scales may require 
more than improving spatial resolution. The use of N- 
body computing in cosmology is possibly overstated as 
is indicated by the differing results obtained when using 
different N-body codes on the scales of interest (Melott 
et al. 1997; Splinter et al. 1998). Furthermore, partic- 
ular features show that caution is in order: (i) the va- 
lidity of the mean field approximation for the gravita- 
tional field strength that is commonly used is question- 
able, especially on the scales of galaxy halos; (ii) the pos- 
sible emergence of soliton states (Gotz 1988) that arise 
as a result of the nonlinear interaction between grav- 
ity and pressure-like forces; solitons have special stabil- 
ity properties and may dominate large scale structure at 
late times; (iii) the behavior of the gravitational field 
strength near high-density regions: from a generic inte- 
gral of the field equations (Buchert 1993a, Sect. 6.1) one 
clearly infers proportionality of the field strength and the 
density, whereas Lagrangian perturbation schemes (in the 
regime where they match N-body runs) show smooth and 
moderately increased field strengths when crossing high- 
density regions. In the extreme case of infinite resolution, 
Zel'dovich's approximation produces caustics in the den- 
sity field, whereas the field strength remains smooth, al- 
though it should blow up at caustics; whether N-body 
codes treat this correctly is questionable in view of the 
strongly varying fields on the spatial as well as temporal 
resolution scales. 

In this situation, analytical models that are capa- 
ble of accessing non-linear scales should be further ex- 
plored. This task is not easy, since — as explained above 
- such models may not be fully guided by compar- 
isons with N-body results. The pioneering model for ad- 
hesive clustering has been built in view of the short- 
comings of Zel'dovich's approximation predicting struc- 
ture decay after their formation: a Laplacian forcing has 
been added ad hoc to the evolution equation for the 
peculiar-velocity. Such a force has the right property 
of holding structures together after their formation and, 
moreover, the model can be solved exactly in terms of 
known solutions of the 3D Burgers' equation (Gurbatov 
et al. 1989; Weinberg & Gunn 1990; Kofman et al. 1990, 
1992, 1994; see also: Gurbatov et al. 1983, 1985, 1991). 
The phenomenon of non-dissipative gravitational turbu- 
lence is mimicked by the so-called "burgulence" (Frisch & 
Bee 2001) and roughhewned into an "effective viscosity". 
However, the performance of the 'adhesion approximation' 
does not dramatically improve the Lagrangian perturba- 



tion schemes, if the latter are subjected to optimization 
techniques as mentioned above: even for a strongly hier- 
archical power spectrum (power law index of n = — 1) 
the optimized Lagrangian schemes yield better cross- 
correlation statistics for the density fields, compared with 
numerical simulations (Buchert 1999). However, compar- 
isons have been conducted only on the level of large-scale 
structure (e.g., Sahni et al. 1994 investigate the evolu- 
tion of voids); the role that the 'adhesion approximation' 
could play on subcluster scales has not been explored. 
Other proposals to analytically model the weakly non- 
linear regime tried to improve the performance for large- 
scale structure, e.g. "freezing" the initial streamlines of 
the fluid keeping the velocity potential constant (so-called 
'Frozen Flow Approximation', Matarrese et al. 1992), or 
treating the gravitational potential as constant ('Frozen 
Potential Approximation', Brainerd et al. 1993; Bagla & 
Padmanabhan 1994). However, these approximations, by 
their nature, are unable to model a realistic evolution be- 
yond the epoch when the time-dependence of the velocity 
or the potential is important, cf. the systematic compar- 
ison of different models carried out by Sathyaprakash et 
al. (1995). 

Lying at the physical core of the problem of cover- 
ing multi-streaming, the 'adhesion approximation' is a 
promising model due to a strong support given to it by 
Buchert & Dommguez (1998): a Laplacian forcing appears 
naturally in a physical description of adhesive gravita- 
tional clustering and the 'adhesion approximation' can be 
derived on the basis of a set of — however, strongly sim- 
plifying — assumptions. The Laplacian forcing does not 
represent a true viscosity, it arises by combining the gravi- 
tational field equations with a dynamical pressure gradient 
due to multi-stream stresses. Consequently, momentum 
is conserved and energy is not lost: the model is time- 
reversible. Multi-streaming implies a reshuffling of energy 
components: bulk kinetic energy is transformed into in- 
ternal kinetic energy, and the gravitational potential en- 
ergy of a high density region gradually dominates over the 
tidal interaction with its environment, thus effectively iso- 
lating the system and — in idealized cases — establishing 
a dynamical equation of state as a relation between gravi- 
tational potential energy and internal kinetic energy. 

This reasoning can be made mathematically precise 
by deriving the equations of motion for the coarse-grained 
(or filtered) density and velocity fields. Given a smoothing 
length L, the system is divided into the degrees of free- 
dom of the scales above and below L, respectively. The 
evolution of the smoothed density and velocity fields is 
in general dynamically coupled to the degrees of freedom 
below L. The standard 'dust model' is recovered by as- 
suming that the dynamical effect of this coupling is neg- 
ligible. When the 'dust model' develops singularities, this 
assumption breaks down, and an improved model has to 
be used that accounts for multi-streaming and the cou- 
pling to the small-scale degrees of freedom. The model 
introduced by Buchert & Dommguez (1998), which we 
call the Euler-Jeans-Newton ('EJN') model, identifies the 
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main source of the corrections to 'dust' as the gravita- 
tional multi-stream (GM) effect. This gives rise to a stress 
tensor of purely kinetic origin — as in an ideal gas - 
which is added to the large-scale gravity already consid- 
ered by the 'dust model'. This stress may be idealized 
in the 'EJN model' phenomenologically as isotropic, and 
modelled with a pressure p{g) depending only on the den- 
sity g, i.e. a dynamical equation of state. 

Numerical tests (Dommguez 2003, Dommguez & 
Melott 2004) show that a polytropic equation of state, 
p oc g 2 ^ 11 in a certain range of densities is characteris- 
tic for a variety of initial conditions. The "anomaly" rj 
measures the deviations from the naive virial prediction 
p oc g 2 , depends itself on initial conditions (77 ranges from 
to 1 as the small-scale power is reduced), and can be 
assigned to the fact that the fluid elements are not iso- 
lated and in a stationary state, as required when deriving 
p oc g 2 from the virial theorem. An exact integral of the 
form p oc g 5 / 3 can be deduced theoretically under certain 
strong assumptions (Buchert & Dommguez 1998). 

Recent developments of the 'EJN model' concern the 
Lagrangian linear regime (Adler & Buchert 1999), where 
solutions can be found by extrapolating known results 
of the Eulerian linear regime; the Lagrangian scheme 
has been developed to second order (Morita & Tatekawa 
2001, Tatekawa et al. 2002) and recently to third or- 
der (Tatekawa 2005); numerical tests employing N-body 
and hydrodynamical simulations have been conducted by 
Tatekawa (2004a,b). 

While the foundations of the 'EJN model' can be 
derived on the basis of the velocity moment hierarchy 
of the commonly employed Vlasov-Poisson system, its 
practical implementation needs (phcnomenological) clo- 
sure conditions to truncate the hierarchy. An attempt to 
systematically go beyond phenomenology is the Small- 
Size Expansion ('SSE') introduced by Dommguez (2000). 
The mode-mode coupling in the equations for the fields 
smoothed over the scale L is estimated under the assump- 
tion that the most important contribution to this coupling 
comes from inter-mode coupling at scales > L. The sys- 
tematic nature of this scheme is due to a formal expansion 
of this coupling in powers of L. The 'dust model' is recov- 
ered as the lowest order term (formally setting L = 0). The 
corrections also yield a Laplacian forcing as well as the 
'adhesion model' under certain simplifying assumptions. 
However, the corrections account not only for the GM 
effect like the 'EJN model', but also for the gravitational 
influence of the density inhomogeneities on scales below L: 
the stress tensor has an additional potential-energy con- 
tribution that can be viewed as a correction to the mean 
field approximation for the gravitational field-strength 
and also provides a Laplacian forcing, although subdomi- 
nant relative to the one due to the GM effect. This effect 
may become dominant if one aims to describe structure 
evolution on galaxy halo scales. The stress tensor in the 
'SSE model' is also a source of vorticity by tidal torques 
and shear stretching. This has been studied recently with 
the Eulerian perturbation expansion (Dommguez 2002). 



In the present work the 'dust', the 'EJN' and the 'SSE' 
models will be presented in a unified framework, as dif- 
ferent closure approximations to a hierarchy of equations, 
and a clearcut description of the weakly non-linear regime 
of adhesive gravitational clustering will be offered. It is ex- 
pected that the gravitational multi-stream process ( GM- 
effect for short) as well as stresses arising from corrections 
to the mean field approximation may help to establish 
the stationary configurations that are usually studied in 
stellar systems theory, where "virialized objects" are, how- 
ever, considered as isolated entities. 

Sect. |2 applies the coarse-graining method to cosmo- 
logical structure formation in the phase space of N parti- 
cles, resulting in a continuum description which features 
the Vlasov dynamics as a subcase. Sect.Oldescribes the re- 
sulting fully non-linear evolution equations for the density 
and peculiar-velocity fields in Eulerian space and summa- 
rizes the assumptions that reduce the general problem to 
the key equation of the 'adhesion approximation'. Sect.0] 
reviews the application of the Eulerian perturbative ex- 
pansion and Section that of the Lagrangian perturba- 
tive expansion. Sect. presents new results beyond the 
perturbative regime. Sect. summarizes the results and 
proposes generalizations of the 'adhesion approximation' 
as well as future prospects. 

Notation : Eulerian coordinates are denoted by x, while 
Eulerian coordinates that are comoving with the Hubble 
flow (i.e. the Lagrangian coordinates of the background 
solution) are denoted by q. In both cases, Lagrangian co- 
ordinates are denoted by X (note that Lagrangian co- 
ordinates q for the homogeneous solution are regarded as 
Eulerian coordinates in an inhomogeneous setting) . Greek 
indices refer to particles, while Latin indices refer to vector 
components; a repeated index indicates summation. 

2. Dynamical equations through coarse-graining 

We consider a system of N identical particles evolving un- 
der gravitational forces. In the cosmological context, these 
particles are the constituents of the dark matter. The pur- 
pose is to model the dynamical clustering of these parti- 
cles, which is supposed to lead to the observed large-scale 
structure. In the Newtonian approximation, the equations 
of motion for the position and velocity of the particles in 
phase space are 

x {a) =v( a \ v {a] =g^ a \ (1) 

where the gravitational field strength is determined 
self-consistently by the Newtonian field equations, 

N 

V («) • = A - 47rGm ^ S(x^ - x { ^) , 

P^a ( 2 ) 

V (q) x = . 

m is the mass of a particle, G denotes Newton's gravita- 
tional constant and A the cosmological constant. 
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There is an alternative way of writing these equations. 
We define the Klimontovich density in the one-particle 
phase space as: 



JV 



f K (x,v, t) := S ( x - x(a) ) 6 (v - v^) . (3) 

a=l 

The equations of motion (^-|2J) lead to a dynamical equa- 
tion for this density 1 : 



df K , df K , df K 



Of 



dx dv 



(4a) 



V ■ g := A - 4nGm J dv f K (x,v,t) , V x g := . (4b) 

This description of the system is exhaustive and much 
more detailed than one is usually interested in. A 
coarser description may be obtained by smoothing the 
Klimontovich density: let £, V denote respectively the spa- 
tial and velocity smoothing scales; then a smooth phase- 
space density can be defined as: 



f(x,v,t) := 



!%^w(^)w(^-)f K {x',v',t) , 



(5) 



where W(-) is a rotationally symmetric coarsening window 
normalized to unity. We can write the evolution equation 
of this new density according to Eqs. (@J: 

<V +v .dl + g.dl = _JL. s (v)_d_ .g (9) ; (6a) 
dt dx dv dx dv 



V. S :=A-4. Gra /*^,„, t ),Vx 9 :=0, (6b, 



with 

S {v \x,v,t) 



£ 3 V 3 I C I \ V 



x («' -v)f K {x',v',t) 



(6c) 



field gravity on the smoothing scales. Notice that Eq. JJJJi) 
is, like the exact Eq. l@Ji), a conservation equation in the 
one-particle phase space, since the total number of parti- 
cles N — J dx dvf(x,v,t) = j dx' dv' Jk(x' ,v' ,t). 

Eqs. jntjb) do not form a closed system of equations 
unless the sources S^', S^ 9 ' can be represented or approx- 
imated as functionals of f(x,v, t). The form and success of 
the approximation depends in general on the initial con- 
ditions and the choice of the coarsening scales C, V. A 
frequently used approximation in many examples of cos- 
mological and astrophysical interest consists of neglecting 
small-scale departures from the coarse variables, « 0, 
~ 0, so that Eqs. JBJi, b) reduce to the Vlasov-Poisson 
system of equations for a smooth f(x,v,t) (e.g., Peebles 
1980; Binney & Tremaine 1987). 

However, to describe the formation of cosmological 
structures, the continuous phase-space density f(x,v,t) 
is still too detailed, since the observables we are interested 
in are the large-scale density and velocity fields. A mass 
density and a mean fluid velocity can be defined by the 
velocity moments of f(x,v,t): 



p(x,i) 



dv f(x,v,t) = ^J2 W 



x — X 

T 



(a) 



(7a) 



,(°0 



c 



pv(x, t) :— m J dv v f(x, v, t) — ^ W 

a 

(7b) 

where the last expressions follow from inserting succes- 
sively the definitions of /, Eq. JSJ, and fx, Eq. ©. 
Eq. i jTb ) defines the Eulerian mean fluid velocity v. 

The evolution equations for these two fields, express- 
ing mass and momentum conservation, follow immediately 
from Eqs. © or from Eqs. 

^+V-(pv)=0, (8a) 

^ + («.V)« = g + i(^-V-7»), (8b) 
dt p 



,(«) 



S «(*, M ):=/$^> 



£ 3 V 3 
x [g{x',t) 



X — X 



W 



V — V 



C J \ V 
g(x,t)]f K (x',v',t). 



(6d) 



The field g(x,t) defined by Eqs. JHt) is the gravita- 
tional mean field associated with the phase-space density 
f(x,v,t). The terms and S^ 9 ' represent the dynam- 
ical coupling to the degrees of freedom removed by the 
smoothing procedure: S^ v ' accounts for the velocity dis- 
persion, while represents departures from the mean 

1 We have disregarded the infinities arising from the evalu- 
ation of the field g(x,t) at the particles' positions. These are 
artifacts of using point particles and are treated by smearing 
the particles over an infinitesimal size. For our formal manip- 
ulations, it suffices with the implicit understanding that the 
infinities can be dropped when they appear. 



V • g := A - AnGp , V x g ■= . (8c) 

g is identical to the gravitational mean field strength de- 
fined in Eq. JSt)- Two new fields have emerged: 



Ti{x,t) 



dv s\ 9 \x,v,t) 



c 



ti a) - §i{x,t)] , (9a) 



i{x,t) := mfdv - Vi(x,t)][vj - Vj(x,t)] f(x,v,t) 
+[v j -v j (x,t)]S$ v) (x,v,t)} 



m \ - 



x — X 

T 



(«) 



r («) (a) 



ViVj(x,t)] . (9b) 
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These new fields, in analogy to and S^ 9 \ represent 
the coupling of the dynamics of the relevant fields p and v 
to the configurational and kinetic degrees of freedom that 
have been averaged out. The symmetric second-rank ten- 
sor field V accounts for the velocity dispersion around the 
mean velocity. As Eq. ^jp) shows, there are two contribu- 
tions when working at the level of the smoothed phase- 
space density f(x,v,t), which are subsumed in a single 
term at the description level of the point particles. The 
vector field J- contains the deviation from mean field grav- 
ity — alternatively, g is the gravitational field created by 
the monopolc moment of the mass distribution in the sub- 
system defined by the window, while is the contribution 
from the higher mass-multipoles generated by the small- 
scale structure. 

The set of hydrodynamic-like Eqs. © can be de- 
rived through the intermediate step of a kinetic equation, 
Eq. or directly from the microscopic particle dynam- 
ics, Eqs. JJJ. In the latter case ,V disappears from the 
definitions J3 El after integrating over the velocities, in- 
dicating that the precise value of V is unimportant for 
Eqs. ©. 

An advantage of the formal derivation through coarse- 
graining as presented here is that it makes explicit the 
smoothing scales C and V, which are usually implicit in 
most applications and not clearly presented. Although 
these two scales are arbitrary, their choice is usually dic- 
tated by the physics of the problem at hand. Typically 
V = 0, and we know of only one application where V ^ 0: 
the numerical codes to solve the Vlasov equation neces- 
sarily have a finite resolution in phase-space and coarse- 
graining in velocity is a method to treat the filamenta- 
tion problem in velocity space induced by this equation 
(Klimas 1987). (In other works (e.g. Lynden-Bell 1967, 
Shu 1978) phase-space coarsening is invoked but the pre- 
cise value of V, being irrelevant, is not addressed). 

The spatial scale £ is usually chosen in concordance 
with the length scale of the phenomena of interest. In the 
regimes of validity of the paradigmatic examples of ki- 
netic theory (Boltzmann equation and Vlasov equation), 
C is much larger than the mean interparticle distance I, so 
that f(x,v,t) does not resemble the spiky Klimontovich 
density fx(x,v,t) (Spohn 1991, and refs. therein). The 
Boltzmann equation for dilute gases with a short range a 
of interaction follows from Eq. (J^Ji) when the mean field 
force g is negligible compared to the effect of the term 
S^ 9 \ which is dominated by binary "close encounters" 
of uncorrelated particles. This holds in the scaling limit 
a — > with I ~ cr 2 / 3 (dilute limit) and C ~ a (continuum 
limit), and the mean free path, A ~ £ 3 /<r 2 , remains finite. 
The Vlasov equation describes the dynamical evolution of 
f(x,v,t) when the interaction is long-ranged and weak. 
For the gravitational interaction, Eqs. (J5J), this holds in 
the scaling limit m — ► with I ~ m 1 / 3 (finite mass den- 
sity) and L ~ m° (continuum limit), so that the particle 
distribution is statistically homogeneous on scales below 
L and the mean field g dominates over the effect of "close 
encounters" contained in S^ 9 \ 



Ma & Bertschinger (2004) study the Klimontovich 
equation for particles moving in a given stochastic large- 
scale gravitational field, intended to be a model of parti- 
cle evolution in a galaxy halo environment; after averaging 
over realizations of the large-scale gravitational field a new 
kinetic equation for the ensemble averaged / is obtained 
that differs from the Vlasov equation (this ensemble av- 
eraging is, in the language of statistical physics, similar 
to an average over 'disorder' rather than an average over 
'thermal noise'). For comparison, our Eqs. (© describe the 
evolution of a single realization of a coarse-grained dis- 
tribution (note that the mathematical structure of both 
averaging approaches are similar with reinterpretation of 
the window function W(-) as a probability density). Of 
course, in addition to coarse-graining one could imple- 
ment ensemble averaging. The addition of a noise term as 
a model of small-scale degrees of freedom involves interest- 
ing physics (e.g., Berera & Fang 1994; Barbero et al. 1997; 
Dommguez et al. 1999; Buchert et al. 1999; Matarrese & 
Mohayaee 2002; Antonov 2004). 



3. Dynamical equations for cosmological structure 
formation 

In this section we apply the method of coarse-graining 
to the case of structure formation in cosmology by 
gravitational instability. In this case it is convenient to 
introduce 'comoving' Eulerian coordinates attached to 
a homogeneous-isotropic solution (Friedmann-Lemaitre 
backgrounds), characterized by the expansion factor a(t) 
(and Hubble's function H = a/ a), with tiny inhomo- 
geneities superimposed as the seeds of structure for- 
mation. In order to subtract the homogeneous-isotropic 
motion, one defines comoving positions q^ a \ peculiar- 
velocities it( Q ) and gravitational peculiar-accelerations 
as follows: 



1 



a 

u {o.) ._ v (a) _ Hx ( a ) ^ 

W W :=g (°) + 47rGgg - A xW 



(10) 



where Qu(t) is the homogeneous background density, sat- 
isfying qjj + 3Hqh = 0. Friedmann's equation 3- = 
A — AuGqh has been used in l|10|) . In terms of these vari- 
ables, Eqs. 111121) are: 



a 



m 

^3 



N 



.08)1 



QH 



V^ Q) x = 



(11) 
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We can now repeat the derivations of SJ21 W e define a 
mass density and a mean peculiar-velocity as follows: 



g(q,t) :-- 



(aLy 



gu(q,t) := 



(aLy 



J2 W 

a. 

J2 W 



(a) 



Q Q 

L 



(12a) 
(12b) 



where L = C/a is the (comoving) coarsening length. This 
length is in principle arbitrary, its use being motivated 
by the features of cosmological structure formation one is 
interested in. In general, one should demand that L will 
be much larger than the (comoving) mean interparticle 
distance, since dark matter discreteness is cosmologically 
irrelevant. In the rest of the paper, when we speak about 
"small/large scales" , we use L to set the comparison scale. 

The equations obeyed by these fields can be derived 
as before. Note that in the following we drop the overbar 
used to distinguish the microscopic velocities and acceler- 
ations from the mean velocities and mean field strengths; 
we always refer to the following equations, so that we no 
longer need this distinction (the time-derivative is taken 
at constant q and L in these equations; see Appendix lAl 
for the relation to Eqs. JSJl): 



dg 
dt 



3Hg 



1 



V q • (gu) = 



(13a) 



^ + -(u-V„)u + Hu = w + - (F--V g -n) , (13b) 
at a g \ a ' 



Vq • w := —4irGa(g — g^) , 



Vq x w 



(13c) 



(13d) 



Here, w (with overbar omitted) is the mean field strength 
of the peculiar-gravity, and we have the "peculiar coun- 
terparts" to the fields ©: 



Fi(q,t) :-- 



Uij(q,t) := 



(aL} 



J2 W 



(a) 



(aL)'- 



J2 W 



Q Q 

L 



q - qW 

z 



(14a) 



(14b) 

If one computes dynamical equations for the new fields 
F and II, one obtains further new fields, and so on. In 
order for Eqs. (|13fl to become useful, one has to close this 
infinite hierarchy by resorting to physically motivated ap- 
proximations that yield the fields F and II as functionals 
of g and u. That a successful closure exists is not guaran- 
teed and this usually depends on the choice of the scale 
L and the initial and boundary conditions. In the context 
of fluids dominated by short-range interactions, there is 
a set of closed hydrodynamic equations when the fields g 



and u vary spatially only over length scales much larger 
than any microscopic length (mean free path, correlation 
length, etc.). In such cases, one must take L to lie be- 
tween these microscopic lengths and the scale of varia- 
tion of the fields. Closure is achieved with the assump- 
tion of local quasi-equilibrium: the coarsening cells of size 
~ L are idealized as systems in almost thermal equilib- 
rium, the corrections to equilibrium being proportional 
to the relatively small spatial gradients of the fields 2 . In 
this manner one obtains the Navier-Stokes and Fourier 
equations. The existence of a hydrodynamic regime where 
these equations hold can be supported theoretically for di- 
lute gases described by the Boltzmann equation by means 
of the Chapman-Enskog expansion (Chapman & Cowling 
1991); for dense fluids, the hydrodynamic regime has the 
status of an experimental fact. 

In the context of cosmological structure formation one 
cannot follow this argument, since the notion of thermal 
equilibrium is not well-defined. We instead discuss three 
closures: 

1. The 'dust model': one assumes that the large-scale 
dynamics is dominated by the large-scale forces and ne- 
glects the coupling to the small scales altogether: one sets 
Fi = 0, Hij.j = in Eqs. I|13fl . an d obtains the Euler- 
Newton system for the peculiar-fields: 



|f +3H e +~\7 q -( Q u) = 
at a 



du 1 . _ . 
-dt + a iU -^ 



u 



Hu = w 



Vq • W 



-A_-nGa(g - g H ) 



Vq X W = 



(15a) 



(15b) 



(15c) 



(15d) 



This system of equations has a long history and has been 
thoroughly studied in the cosmological context (see, e.g., 
Peebles 1980; Sahni & Coles 1995; Ehlers & Buchert 1997; 
Bernardeau et al. 2002). Coming from the context of hy- 
drodynamics, the matter model 'dust' is also called a 
'pressure-less fluid'; from our reasoning it is clear that 
this model is characterized by the assumptions that (i) 
small-scale inhomogeneities are irrelevant, so that mean 
field gravity is dominant, Fi w 0, and (ii) multi-streaming 
is absent and small-scale kinetic degrees of freedom are 
subdominant too, so that "pressure" drops, Hy ps 0. 3 
The 'dust model' was originally applied to the 'top-down 
scenario' of structure formation, e.g., 'Hot Dark Matter 
cosmogony', in which large-scale density inhomogeneities 



2 Notice that this closure requires one to consider the energy 
or, equivalently, the thermodynamic temperature as another 
fundamental field on equal footing with g and u. 

3 According to Eq. 113b ). the same model is recovered if Fi — 
o _1 Ilijj ~ 0, i.e., when the corrections to the mean field are 
counterbalanced by velocity dispersion — we have not explored 
the physical meaning of such a balance condition. 
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grow while the small scales remain homogeneous. The 
model, however, turned out to be a relatively good descrip- 
tion also for the 'bottom-up scenario', e.g., 'Cold Dark 
Matter cosmogony' (Pauls & Melott 1995), in which the 
small scales are always strongly inhomogeneous. Strictly 
speaking, the 'dust model' is not properly denned if there 
is too much initial small-scale power, because of the emer- 
gence of singularities at arbitrarily short times, so that 
a small-scale cutoff is required (this becomes manifest 
in formal perturbative expansions too (Valageas 2002; 
Bcrnardeau et al. 2002)). Our derivation of the 'dust 
model' clearly demonstrates that the smoothing length L 
is indeed a denning ingredient of the model. However, this 
fact has not been properly emphasized in the literature 
- the smoothing length is likely irrelevant in the 'top- 
down scenario', where the initial conditions have a built-in 
length of smoothness (free-flight scale of the dark matter 
particles), but this is not so in a 'bottom-up scenario': in 
this latter case, it was found that a much better agreement 
with N-body simulations is achieved if the initial condi- 
tions are first smoothed, e.g., the "Truncated Zel'dovich 
Approximation" (Coles et al. 1993). 

The 'dust model' correctly describes many features of the 
formation of structure by gravitational instability. It has 
also an important shortcoming, which is the focus of the 
present work, namely that it generates caustics, i.e., den- 
sity singularities, as well as multi-streaming regions (shell- 
crossing), where the velocity field is multi-valued. This 
points to a breakdown of the approximation, so that the 
term Fi — a _1 llyvj in Eq. (|13b ) is no longer negligible. 
In the present work we address two approximations be- 
yond the 'dust model' in which this term is modelled as a 
function of g and u. 

2. The Euler Jeans Newton ('EJN') model: moti- 
vated by hydrodynamics and the theory of stellar systems, 
we have the following phenomenological closure: the cor- 
rections to mean field gravity are neglected, and the veloc- 
ity dispersion is approximated by an isotropic tensor field 
which, in addition, is commonly modelled as a function of 
the local density: 

Fi = , Ily = p{g)5 ij , (16) 

with a given function p(g) > 0, a kinetic pressure. 
(Appendix [U] discusses some examples of the functional 
dependence p(g).) This velocity dispersion can be a conse- 
quence of the multi-streaming or it can be already present 
in the initial conditions. This model is applicable to sit- 
uations in which deviations from the mean field can be 
ignored. For cosmological structure formation in a 'top- 
down evolution', it would be safe to ignore deviations from 
mean field, since then the scales below the collapsing ones 
are relatively homogeneous (e.g., due to the free-flight of 
the dark matter particles as in a 'Hot-Dark-Matter cos- 
mogony'), and the main correction to 'dust' is due to ve- 
locity dispersion. On the other hand, if a scenario that 
entails 'bottom-up' structure formation is investigated re- 
alistically, especially down to galaxy halo scales, then the 



scales below the collapsing ones are very inhomogeneous 
due to the previously formed clusters, and both velocity 
dispersion and stresses due to deviations from mean field 
gravity can be expected to be relevant. 
Given the phenomenological character of the approxima- 
tion (|16f) . one can think of the relationship p(g) as "tem- 
plates" modeling the overall effect of both velocity disper- 
sion and departure from mean field gravity as long as the 
mathematical analysis of Eqs. Ijl7jl is involved. The only 
exact constraint required by this interpretation is that F 
can be written as the divergence of a stress tensor. This 
is the case in the other closure ansatz to be discussed (see 
Eq. JUJl). 

Inserting the ansatz for Fi and into Eq. I|13b ). one 
gets (A q is the Laplacian operator in the variable q): 



8g 
dt 



3Hg+-W q ■ [gu) = 



(17a) 



du 1 . _ , 

— + -(U- V q )u + Hu 

at a 



w V q g 



w 



no 
A 2 (g) 



^-AqW , 



Vg • w = —ATtGa(g — qh) 



Vg X W 



(17b) 
(17c) 
(17d) 



where we have denned a density-dependent length (see 
Eq. I|34[) for the relation to Jeans' length): 



p'(q) 



inGg 



1/2 



(18) 



Eqs. I|17b ,.b) are Euler's equations for a compressible fluid, 
Eq. I|17b) is known as Jeans' equation in stellar sys- 
tems theory. This model was proposed by Buchert & 
Dommguez (1998) starting from Vlasov's equation for the 
one-particle distribution function, so that corrections to 
mean field gravity were neglected from the outset and only 
velocity dispersion remained. In that work it was also sug- 
gested that the phenomenological ansatz Eq. (|16f) could 
be replaced by a different closure condition leading to 
p oc p 5 / 3 (see Appendix [UJ|. 

The case of isotropic stresses formally covers the hydro- 
dynamics of a perfect fluid or, as common in studies of 
stellar systems, polytropic models. However, there is no 
compelling reason for isotropic stresses in collisionless sys- 
tems before the onset of virialization. Even "virialized" 
systems will in general maintain an anisotropic compo- 
nent. At the moment, we consider the isotropy assump- 
tion as a good working hypothesis to understand the role 
of the interaction between multi-stream forces and self- 
gravity: the velocity dispersion ellipsoid is approximated 
by a sphere. (Interesting considerations of the influence of 
the anisotropic part have been reported by Maartens et 
al. 1999; see also Barrow & Maartens 1999.) 
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3. The Small-Size Expansion ('SSE'): this is a 
method proposed by Dommguez (2000, 2002) to formalize 
the notion that the coupling to the small scales is in some 
sense weak. One argues that the corrections Fj and Ilyj 
are determined mainly by the largest scales contributing 
to them, i.e., those close to L. For instance, in the bottom- 
up scenario, the large-scale dynamics is sensitive mainly 
to the motion of the most recently formed clusters as a 
whole, and not so much to their internal structural details 
due to the trapped particles, so that one takes L ~ typ- 
ical size of these clusters (which play the role of effective 
particles of size ~ L). The mathematical implementation 
of this idea leads to a formal expansion of Ft and IL,j in 
powers of the smoothing length 4 L (see Appendix \B\ for 
an outline of the derivation): 



F t ^ BL 2 (d k g)(d k w t )+0(L 4 ) 



lit 



BL 2 Q{d kUi ){d k Uj)+0{L A ) , 



(19a) 
(19b) 



where B > is a dimensionless constant of order unity 
depending on the shape of the smoothing window W(-), 
Eq. I|B.5J) . The correction to mean field, being proportional 
to dkWi, is approximated by the tidal forces induced by the 
mean field gravity. With Eqs. Ijl3h .d). it can be written 
also as the divergence of a second-rank tensor: 



Fi -> BtfdjigdjWi + 2ixGag 2 b %i ) + 0(L 4 



(20) 



The velocity dispersion tensor is represented through the 
kincmatical parts of the peculiar-velocity gradient: expan- 
sion rate 9 := (l/a)V q • u, vorticity u> := (l/a)V q x u, 
and shear av 



IL, -> a 2 BL 2 g 



I- _ 1 



ax 1 

UOii + (Ta + -EijkUJk 



(21) 



1 2 2 

—6 Sij + -9<Jij + (Jik<Jkj — 
1, . 



-{iJdij -LJiUj) 



0(L 4 ). 



(22) 



In the context of the 'SSE model', we recover the 'dust 
model' to order O(L ). To order 0(L 2 ), we obtain the 
following equations: 



du 1 

dt a 
BL 2 



^+3Hg+-V q -(gu) = 
at a 



(U ■ Vq)ti + HU = W 



(23a) 



(23b) 



( (Vqg ■ V q )w - -Vq ■ [g(d k u)(d k u) 



Vq • w = —4:nGa(g — qh) 



(23c) 

4 The denomination "large-scale expansion" used in 
Dommguez (2000) may lead to confusion. 



Vq x w = . (23d) 

The expansions (|19|l are formally exact. The assumption 
underlying the 'SSE' is that the first few terms provide an 
accurate description of the dynamical effect of the mode- 
mode coupling for the problem at hand. For comparison, 
this is not true in the hydrodynamic regime of fluids dom- 
inated by short-ranged interactions, where the high-order 
terms in the expansions have to add up to yield the usual 
L-independent hydrostatic pressure and viscous force. 
The closure assumption (|19b ~l for the velocity dispersion 
has been also used to model the influence of subresolution 
degrees of freedom of the velocity field in "Large-Eddy 
Simulations" of turbulent flow (Pope 2000). It has been 
termed the 'gradient model' (see, e.g., Vreman et al. 1997) 
and it was introduced as part of the 'Clark model' (Clark 
et al. 1979) (a model which includes an extra additive term 
in the expression for Ily ) . 

Both the 'EJN' and the 'SSE' models provide an ex- 
tension of the 'dust model' with the potential to improve 
on it by preventing the formation of singularities. The 
'EJN model' relies on a phenomcnological assumption con- 
cerning the kinetic stress Ily , but has the advantage that 
the resulting equations are well established. In favor of 
the 'SSE' model is the fact that the corrections to 'dust' 
can be obtained systematically beyond phenomenology, 
but the mathematical and physical status of the resulting 
equations is little explored yet. 



3.1. Where do these equations hide the 'adhesion 
approximation ' ? 

In this subsection we get a first picture of the effect of the 
corrections to the 'dust model'. In the 'adhesion model' 
(Gurbatov et al. 1989), the evolution of 'dust' is modelled 
with Zel'dovich's approximation (Zel'dovich 1970, 1973), 
a powerful and mathematically manageable description of 
the exact evolution (first-order Lagrangian perturbation 
solution, see H5.1|l . Accordingly, one assumes in Eq. H15h) 
the proportionality of peculiar-velocity and -acceleration 
fields (c/. Peebles 1980; Buchert 1989, 1992; Bildhauer & 
Buchert 1991; Kofman 1991; Buchert 1993b; Vergassola et 
al. 1994): 



w = F(t)u, F(t) = 4irGg H ^ 



bit) 



b(t) 



(24) 



where b(t) is identical to the growing-mode solution of 
the Eulerian linear theory of gravitational instability for 
'dust' (i.e., it solves the equation b + 2Hb~ AnGgnb = 0). 
Notice that this "slaving" assumption implies an irrota- 
tional flow: V q x u = 0. Changing the temporal variable 
from t to b and defining a rescaled peculiar-velocity field 
ii = u/ab, Eq. I |15b 1 becomes, after inserting Eq. I|24|) . 



du 

Tb= Q 



db 



d_ 

~db 



— := — +«• V ( 



(25) 
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The evolution according to this equation leads (kinemat- 
ically) to the formation of singularities, where g — ► +00, 
and the velocity field becomes multi-valued. Following the 
idea of the 'adhesion model', the corrections to 'dust' are 
assumed to be negligible almost everywhere, i.e., except 
at the localized singularities. The corrections are then es- 
timated with the 'dust' solution (i.e., with Zel'dovich's 
approximation) : 

- The 'EJN model': inserting the parallelism condi- 
tion (|24|l into Eq. Q17bh and performing the same change 
of variables as before, we obtain, instead of Eq. Q25J1 . the 
following corrected equation: 

^ = /lAgii , (V q x u = 0) , (26a) 

with the GM or gravitational multi-stream coefficient 
(Buchert & Dommguez 1998, Buchert et al. 1999) 



M := -^A 2 (g) > , 
a z b 

formally playing the role of a "viscosity" . 



(26b) 



- The 'SSE model': the parallelism condition allows 
us to express w and g (via Eq. J23t)) in terms of u. The 
correction term is hypothesized to be relevant only near 
places where Zel'dovich's approximation yields singulari- 
ties, where |Vu| — * 00, so that it can be approximated by 
neglecting subdominant powers of Vit: 



F 



in, 

a ■ 



BL 1 



-djKVq ■ u)(d k Ui)(d k Uj)} 



(27) 

To simplify this expression further, one recalls that, 
within Zel'dovich's approximation, density singularities 
arise generically due to a locally plane-parallel collapse 
and the velocity gradient looks locally like diUj w (V q ■ 
u)niTij, where the unit vector n defines the plane of lo- 
cal collapse. In order to compute the correction term, we 
assume that the spatial dependence of n is negligible com- 
pared to the (large) gradient of the fields along n, i.e., as 
if it were a globally plane-parallel collapse, and we obtain: 



1 



-n, 



(V q ■ u)A q Ui 



Therefore, Eq. (j2St>) simplifies, after the change of vari- 
ables, to the same Eq. tyZbh ) with a different GM coeffi- 
cient: 

^.= ^g>0. (28) 
bg H 

The dominant contribution in Eq. (|27|l is the one by ve- 
locity dispersion. It can be easily checked that the correc- 
tion to the mean field gives, under the same assumptions, 
a subdominant but also positive contribution to the GM 
coefficient. 

Eq. Ij26(l resembles the well-known key equation of the 
original 'adhesion approximation' (Gurbatov et al. 1989), 



to which it reduces when /i=constant — > + . For con- 
stant fi it can be solved analytically for curl-free flows 
(Hopf-Cole transformation of the 3D Burgers' equation); 
this solution shows that the singularities predicted by 
Zel'dovich's approximation, Eq. (f2~5|l . are indeed regular- 
ized by the /i-term (see, e.g., Vergassola et al. 1994); 
for an interesting application of Burgers' equation and 
further insight see also the model of Jones (1996) for a 
two-component system. In the more general cases of a 
density-dependent GM coefficient, no analytical solution 
was found. Nevertheless, the application of boundary- 
layer theory (Dommguez 2000) shows that density sin- 
gularities are still regularized (provided p(g — > +00) ~ g 1 
with 7 > 1 in the case of the 'EJN model', so that pres- 
sure can oppose gravity successfully). The solutions for 
the velocity and density fields are qualitatively identical 
to those of the original 'adhesion model', developing into 
a shock structure of infinite density in the limit fi — > + . 

This derivation of the 'adhesion approximation' offers 
insight into the physics involved. The /x-term, which was 
phcnomenologically motivated by analogy with viscosity, 
is not related to a truly dissipative process. From a mathe- 
matical point of view, the starting point ('EJN model' l|17fl 
or 'SSE model' Q23J1). the approximations, and the final 
Eq. (|26[) are formally time-reversible, i.e., invariant under 
the transformation t — > —t, u — > —it. The correction to 
'dust' can transform mean kinetic energy (oc £>|m| 2 ) into 
internal kinetic energy (oc ffj and conversely (also into 
internal gravitational potential energy via Fi, but, as we 
have shown, the 'SSE model' predicts this contribution to 
be subdominant near density singularities in the present 
approximation). The tendency to compression in collaps- 
ing regions favors the correction to behave as a drain of 
mean kinetic energy, mimicking viscous dissipation; it is 
easy to check that the correction supplies energy in ex- 
panding regions (see H6.4[l . whose expansion is thus ac- 
celerated. This mechanism seems to be universal: this is 
the reason that "adhesiveness" arises in the two differ- 
ent models we study, 'EJN' and 'SSE', and that it can be 
modelled in a qualitatively correct manner by the simpli- 
fied 'adhesion model' /i = constant. The GM coefficient 
depends on time and on density, and this leads to differ- 
ences between models concerning the inner structure of 
the regularized density singularities. A density indepen- 
dent GM coefficient (as in the original 'adhesion approxi- 
mation') is obtained with the imposed relationship p oc g 2 
corresponding to the naive application of the virial condi- 
tion, while the 'EJN' and the 'SSE' models coincide when 
p oc g 3 . 

3.2. Relaxing key assumptions 

Summarizing, to derive Eq. (|26|l one makes two important 
assumptions: 

1. The mean peculiar-gravitational field w for the evolu- 
tion of 'dust' is modelled by a first-order Lagrangian 
perturbation solution which includes parallelism of 
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peculiar-velocity and -acceleration (Zel'dovich's ap- 
proximation). 

2. The so-approximated evolution of 'dust' is used to 
estimate the correction to 'dust', which for consis- 
tency must be assumed to be negligible everywhere ex- 
cept where singularities would develop. (For the 'SSE 
model', one makes explicit use of the additional con- 
straint of a locally plane parallel singularity, which 
is nevertheless a consequence of Zeldovich's approxi- 
mation and hence not expected to introduce new fea- 
tures.) 

Relaxing (1) means using a better approximation to the 
exact evolution of 'dust'. Zel'dovich's approximation, orig- 
inally designed as an extrapolation of the standard linear 
theory of gravitational instability into the weakly non- 
linear regime, features "slaving" of peculiar-velocity and 
-acceleration. This simple, but physically well-motived 
condition, if assumed initially, is maintained up to sec- 
ond order in Lagrangian perturbation theory (see: Buchert 
1989, 1992 for the general first-order solution, Buchert 
& Ehlers 1993 for explicit discussions within the gen- 
eral second-order solution, and Buchert 1994 for the cor- 
responding third-order solution). Although corrections 
to Zel'dovich's approximation in general break the par- 
allelism of w and it, the departures are not expected 
to be large in the perturbative regime, given that the 
'Zel'dovich approximation' performs well before shell- 
crossing. Nevertheless, using better approximations for the 
evolution of 'dust' will not prevent singularities from ap- 
pearing, because this is a general property of the 'dust 
model' (see UHJ. 

While relaxing (1) implies construction of sophisti- 
cated non-perturbative approximations (since only then 
can we expect to substantially improve on the perfor- 
mance of Lagrangian perturbation schemes, optimized to 
match N-body simulations for any kind of dark matter), 
assumption (2) appears more like an unnecessary simplifi- 
cation. It is crucial for the derivation of equations similar 
to the key equation of the standard 'adhesion approxima- 
tion', but it must be relaxed if one intends to improve 
over the 'adhesion model' in order to describe the sub- 
sequent dynamical evolution inside the collapsing struc- 
tures. Assumption (2) may be rephrased by saying that 
the evolution is dominated by convection, (u ■ Vq)M, over 
"viscosity" , \xA q u — in the hydrodynamic jargon, this 
is the limit 'Reynolds number — > +oo'. Relaxing this as- 
sumption, i.e. finite Reynolds number, means taking ac- 
count of the back-reaction of the correction to 'dust' on 
the trajectories of fluid elements away from singularities 
too. In general, this will also imply corrections to the par- 
allelism of w and u as well as those due to the exact 'dust' 
evolution (see, e.g., Eq. j55)l). 

Menci (2002) has studied the correction to the origi- 
nal 'adhesion approximation' when assumption (1) is re- 
laxed. The starting point is the set of Eqs. I|13[) with the 
correction to 'dust' already modelled as fiA q u, \x = con- 
stant — > 0. The correction to parallelism, w — Fu, is es- 



timated perturbatively in the limit of small velocities or 
times (the initial condition is taken to satisfy parallelism) 
— the procedure employed by Menci is not an expansion 
in the "viscosity" /x, although it might formally appear so. 
The results agree somewhat better with N-body simula- 
tions concerning the details of the forming structures at 
small scales. 

The solution to the 'adhesion model', Eq. (|26h ) with 
constant /x, cannot be computed as a regular expansion in 
/i, because the naive zeroth-order term, given by Eq. (|25(l . 
is undefined after shell-crossing. The exact analytical so- 
lution to this model demonstrates that the limit /i — ► 
yields a weak solution of Eq. I(25|l . exhibiting discontinu- 
ities in u. In the pedagogical review by Frisch & Bee (2001) 
it is illustrated how much the Lagrangian solution of a 
'multi-streamed dust model' differs from the Lagrangian 
solution to the 'adhesion model' after shell-crossing. In 
particular, it is explained how the emergence of multi- 
streaming cannot be captured by a Taylor expansion in 
the time variable of the solution to Eq. 125(1 . As a matter 
of fact, the corresponding perturbative expansion of the 
solution to the Vlasov-Poisson system is identical to the 
one of the 'dust model' when the initial condition is ex- 
actly single- streamed (Valageas 2001) 5 . This problem can 
be overridden if a small amount of multi-streaming is al- 
lowed in the initial conditions, e.g. in the form of velocity 
dispersion. 

Ribeiro & Peixoto de Faria (2005) also attempted to 
provide a physical motivation for the 'adhesion model': 
starting from the velocity potential tp (with u = \7 q (p/a), 
it is assumed that the complex variable tp — ^/pexp (iip/v) 
obeys a Schrodinger equation, where v is a constant with 
the appropriate dimension. From this assumption, it is 
found that u satisfies Eq. (|13b^l with a correction to 'dust' 
that can be written as a functional of the density. It is not 
obvious that with this correction term one can recover the 
'adhesion model': it cannot be brought into a simple form 
oc A q u after using the parallelism approximation, and the 
reasoning offered by the authors is wrong from the outset 
because their Eq. (31) is algebraically false. 

In the following sections we study the 'EJN' and the 
'SSE' models beyond assumptions (1, 2). First, we con- 
sider the standard Eulerian and Lagrangian perturbative 
expansions, both closely related to the expansion in time 
just mentioned. One of the key points in the derivation of 
the 'adhesion model' is the parallelism relation between 
w and u, Eq. I|24jl. This relation holds at the (Eulerian 
as well as Lagrangian) linear order of the 'dust model'. 
Application of the perturbative techniques to the 'EJN' 
and 'SSE' models will show how this relation is modi- 
fied by the correction to 'dust'. Later we consider non- 
perturbative approaches: this is a mathematically difficult 



5 Bharadwaj (1996) extended this claim to the case of multi- 
streamed initial conditions in the simplified model of free par- 
ticles. This conclusion is, however, false because of a mathe- 
matical mistake in his Eq. (34). 
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task, and we are only able to collect some results and pro- 
vide some hints for future work. 



Then, Eqs. (|31[) can be simplified to a single equation for 
the density contrast 5 := q/qh — 1: 



4. The Eulerian perturbative expansion 

The expansion parameter of the Eulerian perturbative ex- 
pansion is the amplitude of the departures of the initial 
conditions from homogeneity. Formally, one writes: 



Q : 



OO 



(29) 



where e -C 1 is an expansion parameter, which can be 
taken to be the variance of the initial density fluctuations, 
e = er(io), with (note that a is finite because the density 
field has been smoothed over the scale L, Eq. I^Ji)) 



([g(M - QH{t)f 



(30) 



Here, (• • • ) denotes ensemble average over the initial con- 
ditions. Inserting these expansions into Eqs. (|13|l and 
equating powers of e, one obtains a set of linear equa- 
tions for g( n \q,t), u^{q,t). The zeroth order is the ho- 
mogeneous background, Q^(q,t) = gnit), u^°\q,t) = 0. 
Due to the gravitational instability, the fields g( n \ 
(n > 1) eventually grow in time and the expansions l|29Jl 
(likely to be understood as asymptotic expansions) cease 
to be a useful representation of the solution after some 
time T, which can be roughly estimated by the condition 
a{T) = 1. 

4.1. The linear regime, 0(e) 

To this order, one obtains the following equations: 

d 8 V> 



dt 



3#£> (1) + -QH^q ■ U (1) = 



(31a) 



02$(D saw 

+ 2H- 



dt 2 



dt 



AttGqhS^ 



(34) 

where we have used the length defined in Eq. (|18fl . Notice 
that in the absence of cosmological expansion, the time- 
dependence of A induced by qh drops and it reduces to 
the usual definition of Jeans' length (up to a factor 27r, 
e.g. Peebles 1980). Eq. (|3~4*|l can be solved analytically 
for particular choices of the cosmological background and 
the dependence p(g) (Haubold et al. 1991). Qualitatively, 
one finds that the pressure term is irrelevant for inhomo- 
geneities on scales 3> A, which grow according to the 'dust 
model', while those on scales <C A are damped by pressure 
and may exhibit sound-like oscillations. 

The 'EJN model' provides a correction to the par- 
allelism relation already at the linear order. Following 
Buchert et al. (1999) 6 , the solution to Eq. ijSH) and the 
departure from parallelism can be computed as a series in 
powers of A (unlike in the fully non linear problem, the 
limit A — > in Eq. (|34|l is regular). Asymptotically in time, 
we may write (C(t) is a dimensionless time-dependent co- 
efficient): 



i£>« = F(t) 



u 



(1) - ^\ 2 (9H)A q u^ , fti 

n z 



C(A 4 ) 



F(t)uM +c(t)^±V qg ^ + 0(X 4 ) 



CLQH 



(35) 



where the last equality follows by using Eqs. i|31bll8|) . 

The correction is proportional to the extra acceleration 
due to the pressure gradient. Since C(i) ^ 1 in general, 
the peculiar-velocity is not parallel to the total peculiar- 
acceleration (gravity-l- pressure) . 



^ + Hu^ = W P> + J_ ( fW - -v q ■ nw 

dt qh V a 

(31b) 
(31c) 

(31d) 

For the 'dust' and the 'SSE' models, F (1) = and 
n*' 1 -' = 0, and they do not differ at this order. One can 
show that the vorticity, u) = x u, decays with time, 
and that the peculiar-velocity becomes parallel to the 
peculiar-acceleration asymptotically in time, 



Vq • w {1) = -47rGa(p (1) - g H ) 



V„ x 



F{t)u 



(i) 



For the 'EJN model', = 0, but: 
nV=p'(g H )gM5 i: 



(32) 



(33) 



4.2. The non-linear orders, 0(e n ), n > 1 

The Eulerian perturbative expansion beyond the linear 
order has been intensively studied for the 'dust model' in 
recent years (see, e.g., Bernardeau et al. 2002 and refer- 
ences therein). The predictions of the 'SSE model' differ 
from those of the 'dust model' at non-linear orders. An 
important difference concerns the vorticity. With full gen- 
erality (i.e., to all orders of the perturbative expansion), 
one can easily show from Eq. I)13b) by application of the 
vector identity (u-V q )u = (l/2)V q |u| 2 — tix (V,xu) that 
(here, a comma denotes partial derivative with respect to 



6 In that work, only the particular case of a polytropic depen- 
dence, p(g) oc g 1 , in an Einstein-de Sitter cosmological back- 
ground was considered, but it can be checked that Eq. I|M5^ 
holds generally. Also in that work, a minor mistake was made 
and the coefficient C(t) quoted there is missing an extra time- 
dependent factor when 7 > 13/6. 
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comoving Eulerian coordinates) 
duj 1 

+ 2HuJi + - {UjUJtj — Ui^Ldj + UjjLUi) 



g a 



(36) 



For the 'dust model' (vanishing r.h.s.), this equation can 
be integrated exactly along the flow lines, Eq. l(3"§j) (e.g., 
Serrin 1959, Buchert 1992, Eq. (26) ff): 



= (£) 



2 CJj(X,t ) 



dPj(X,t) 
dX 3 



det (da + 



dP,(X,t) 



(37) 



We see that only the corrections to 'dust' can be a source 
of vorticity: if the initial vorticity vanishes (as is usu- 
ally assumed in cosmological models of structure forma- 
tion by virtue of the decay of the linear vorticity), then 
the evolution according to the 'dust model' predicts that 
the vorticity remains zero (Kelvin's circulation theorem). 
This remains true for the 'EJN model', because the pres- 
sure is represented as a function of density only and the 
r.h.s. of Eq. I|36|l vanishes. The 'SSE model', however, 
does incorporate the non-linear generation of vorticity: 
when Eq. (|36|l is evaluated to order 0(e 2 ) with the expan- 
sions l|2"9"|) . one obtains: 



(2) 



Of 



2Huf- ] = BL 2 e llk 



o {1) w {1) 



1 



(U 



(1) „.(1) 



jm 



(38) 

The source terms of vorticity are the tidal and shear 
stretching and can be evaluated as a functional solely of 
the linear density field g^ by resorting to Eqs. l|3"Tb.l3"2"l). 
A detailed study of the growth of the ensemble average 
(w 2 ) and its dependence on the initial conditions can be 
found in (Dommguez 2002). As soon as the velocity field 
has a non- vanishing vortical component, parallelism with 
the gravitational peculiar-acceleration is rigorously ruled 
out. (But this does not mean that deviations from paral- 
lelism must be necessarily large or important). Vorticity 
acts against collapse (see © an d is likely to play a role in 
shaping the internal structure of collapsed regions, which 
in the 'adhesion model' are structureless. 



5. The Lagrangian perturbative expansion 

To perform the Lagrangian perturbative expansion, 
we have to transform the equations with respect to 
Lagrangian coordinates and follow the field quantities 
along comoving trajectories of the fluid elements: 



q = X + P{X,t), u = aP, 



P(X,to) = 0, (39) 



with the Lagrangian time-derivative operator ' := J||x = 
■§l\ q + ~^ ~$q-' (Note that X + P are integral curves of the 
mean peculiar-velocity, scaled by a(t)). X are Lagrangian 
coordinates and P(X,t) is the displacement of the fluid 
volume elements away from their initial positions; this is 



the expansion parameter in the Lagrangian perturbative 
approach. This method has been applied successfully to 
study the 'dust model', Zel'dovich's approximation being 
a subcase of the first-order solution, as well as higher- 
order solutions. It has been recently applied to the 'EJN 
model' too (see the references below). Concerning the 'SSE 
model', there is no difference with the 'dust model' to 
linear order. 

In view of later considerations, for the application of 
the Lagrangian methods we take an unconventional route 
and first rewrite Eqs. (|13|l as follows. We obtain an evo- 
lution equation for the peculiar-gravitational field w by 
eliminating the density in Eq. j!3h) with the help of 
Eq. (fT5h ). By formally integrating the divergence one gets 
(Buchert 



w + 2Hw - AttGqhu = TZ , 
TZ := ~[(u ■ V q )w - u(X7 q -w)+X7 q xT] 



(40a) 



The vector T is essentially fixed by the other equation 
satisfied by w, Eq. ifHfl Ti: acting with V q x on Eq. l|l0"k1 
and subjecting T to the Coulomb gauge condition, X7 q ■ 
T = 0, one finds 

A q T = AnGaVq x (gu) = 4:TrGa[gV q x u + V q g x u] 

= (AirGagn - V q • w) V q x u + u x A q w . (40b) 

That is, T is the vector potential (up to an unimpor- 
tant factor) of the peculiar-current density j := gu (i.e., 
4-7rGaj = — X7 q x T — V q -dt(a<f>), with denoting the 
scalar peculiar-gravitational potential). (If w =: -^V g ^ 
is initially irrotational, Eqs. I|40b.bl propagate this condi- 
tion in time.) Finally, Eq. I|13b ) gives the evolution of it: 



it + Hu 



-v q -n 



(40c) 



Thus, we have replaced Eqs. i|13fl for g, u, and w by the 
new set of Eqs. (|40() for u and w. They hold independently 
of what the correction to 'dust' looks like. These equations 
still retain the nonlocal nature imposed by gravity through 
the vector potential T given by Eq. l|4T)h): we obtain "lo- 
cal" evolution equations only if the source of this equation 
vanishes. Physically this means that the peculiar-current 
density has to be irrotational, and, for irrotational flows, 
that the mean flow follows the gradient of the density 
field 7 . In the following, we will not be further concerned 
with the residual vector field TZ. We just acknowledge that 
it is a functional of the dynamical fields u and w so that 
the system of equations l|40h,-c) is closed for given F and 

n. 



7 If such conditions are imposed on the problem, T becomes 
a harmonic vector function which can be set to zero for periodic 
boundary conditions of the cosmological fields on some large 
scale (since then the only harmonic functions are spatially con- 
stant and can be set to zero due to the invariance of the basic 
equations with respect to spatially constant translations). 
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5.1. The single-streamed case of 'dust matter' 

For the 'dust model' the system of equations l|4()h -c) has 
been studied thoroughly in (Buchert 1989). There, the 
currently used notion of the weakly non-linear regime 
has been defined as "Lagrangian linearization". The co- 
moving trajectory field P of Zel'dovich's approximation 
(Zel'dovich 1970, 1973) can be understood as a subclass 
of the linear solutions in Lagrangian space (Buchert 1989, 
1992). For an Einstein-de Sitter cosmology with initial 
peculiar- velocity u(X, to), the subclass corresponding to 
Zel'dovich's model is P(X,t) = |[&(t) - b(t )]u(X ,t ) t , 
where b(t) — a(t) here. For all background cosmolo- 
gies including "curvature" and a cosmological constant 
Zel'dovich's approximation is given in (Bildhauer et al. 
1992; for the flat models with cosmological constant, see 
also the supplement by Chernin et al. 2003). 

This well-known result can be easily obtained from 
the general system of equations (^Ofa -cl by the following 
reasoning: in the Lagrangian picture the Eulerian dynam- 
ical variables u(q,t) and w(q,t) are replaced by the sin- 
gle dynamical variable P(X,t). Lagrangian linearization 
means that the equations have to be linearized with re- 
spect to P. Eqs. (|39(l and l|4()h) (without corrections to 
'dust') express the fields u and w as linear functions of 
P. It follows that Eqs. I|40b .b) have to be linearized in the 
fields u and w. For simplicity we restrict the initial con- 
ditions to irrotational flow, so that the velocity remains 
irrotational through the 'dust' evolution, Eq. Q37fl. Then, 
Eq. I0TJ)) implies, upon linearization, T — 0. From l|40h,) 
we can therefore immediately drop the manifestly non- 
linear terms in the residual vector field TL. Introducing 
the longitudinal and transversal parts of u and w with 
respect to Lagrangian coordinates, e.g., Vo • w T = 0, 
Vo x w L = 0, where Vo denotes the nabla operator with 
respect to Lagrangian coordinates, the remaining equa- 
tions for 'dust' read: 

w L + 2Hw L - 4ttGq h u l = , u L + Hu L - w L = , 
u T = w T = . (41) 

These equations are satisfied if the inhomogeneous dis- 
placement field P, split into its longitudinal and transver- 
sal parts, obeys: 

P L + 2HP L - AnGg H P L = , (42) 
P T = , 

which can be easily demonstrated by reinserting the defi- 
nitions of the fields u and w in favor of P — P L + P T . If 
we restrict the general solution of Eq. 1)42(1 to the initial 
data w(X,to) = (const./ to) u(X, to) (const. = 1 for an 
Einstein-de Sitter cosmology), we obtain Zel'dovich's ap- 
proximation. The rigorous calculation, starting from the 
fully transformed system of equations, may be found in 
(Buchert 1989, 1992). 



5.2. The multi-streamed case modelled by dynamical 
stresses 

In the regime when the 'dust model' predicts multi- 
streaming, the corrections to 'dust' become relevant. We 
look for the Lagrangian linear form by also dropping the 
residual vector field Tt, (since initial irrotationality is pre- 
served to linear order by the evolution in both the 'EJN' 
and the 'SSE model'). The 'SSE' correction to 'dust' is also 
non-linear, Eq. (|23b L so that this model does not differ 
from the 'dust model' to linear order. For the 'EJN model', 
Eq. l)17bL however, we have to linearize additionally the 
(Eulerian) Laplacian together with the GM coefficient 
along the comoving trajectory field. Using transformation 
tools developed by Adler & Buchert (1999, Appendix A) 
we find (for isotropic II): 

^M AqW = ^l AoW + (P"). (43) 

The complete set of Lagrangian linear equations may then 
be cast into a form that reproduces the 'dust' result l|41|l 
on their left-hand-sides: 

w L + 2Hw L - 4irGg H u L = , 
u L + Hu L -w L = \\ 2 (q h )A w l , (44) 
u T = w T = . 

Again, the rigorous way of deriving the equation obeyed 
by Lagrangian linear displacements is the following: we 
first transform the basic equations into Lagrangian form, 
then linearize the Lagrangian system with respect to P. 
This has been done in (Adler & Buchert 1999) with the 
result that P obeys: 

P L + 2HP L - AkGq h P L = 4nGg H ^T Ll AoP L , 

P T + 2HP T = ; here : P T = . (45) 

After some (in this case more involved) manipulations one 
can find the equations for u and w above by inserting their 
definitions in terms of P into Eqs. (|45fl . The differential 
operator acting on P L is formally identical to Eq. (|34|) . 
and Eq. (|45|l can thus be solved by extrapolation of the 
known Eulerian linear solutions (Adler & Buchert 1999). 

Evidently, the adhesive term consists of a Lagrangian 
Laplacian and, therefore, it involves, in addition to the 
convective non-linearities hidden in the overdot, non- 
linearities in Eulerian space when mapping it back using 
the inverse solution of the mapping q = X + P(X, t). 

In general, parallelism of peculiar-velocity and 
-acceleration is violated. By computing the time deriva- 
tive of the difference w — F(t)u with Eqs. I|45|l . one can 
show that in general it does not decay to zero asymp- 
totically for large times. The correction is proportional to 
Ao« to lowest order in A 2 . This correction is similar to the 
one found in the Eulerian linear approximation, Eq. I(35|) . 
but in terms of the Lagrangian gradients. 
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5.3. Beyond first order 

The Lagrangian perturbation scheme including pressure 
(Adler & Buchert 1999) has been pushed to second or- 
der (Morita & Tatekawa 2001; Tatekawa et al. 2002) and 
recently to third order (Tatekawa 2005). Results on cross- 
correlation statistics with N-body simulations (Tatekawa 
2004a) and a comparison of the corresponding density 
fluctuations with results of hydrodynamical simulations 
(Tatekawa 2004b) has demonstrated that multi-stream 
forces tend to be underestimated by Lagrangian perturba- 
tive models. In the next section we shall investigate why 
this happens. 



6. The non-perturbative regime 

6.1. Exact evolution equations I: general case 

Following a systematic perturbative approach we could 
provide models for adhesive gravitational clustering in the 
weakly non-linear regime. This has allowed us also to re- 
lax the parallelism assumption (|24|l that was a necessary 
ingredient of the standard 'adhesion model' but, as we 
have shown, is not expected to hold in the multi-streamed 
regime. 

Finding an extension into the non-perturbative (both 
Eulcrian and Lagrangian) regimes is an involved math- 
ematical task. Notwithstanding, it should be attempted 
in view of the fact that even the Lagrangian perturba- 
tion approach falls short of capturing the action of multi- 
stream forces, since a genuine property of adhesive mod- 
els in the simplest cases involves an Eulerian Laplacian. 
Thus, we expect the best approximations to be of hy- 
brid Lagrangian/Eulerian type like the standard 'adhesion 
model', meaning a non-linear equation both in Eulerian 
space (due to the convective non-linearities) as well as in 
Lagrangian space (due to the Eulerian gradients in the 
forcing) . 

An exact equation for the comoving deviations of the 
trajectory field, P(X,t) defined by Eq. (|39fl . can be ob- 
tained from Eqs. 140|) : formal integration of Eq. iAOh ) 
yields 



(46) 

with the integration constant chosen such that the 
Lagrangian coordinates are defined by the initial datum 
for w, s ATrGa QH{to)P{X,t ) = w(X,t ). Then, from 
Eq. H4Ub) one obtains the general result: 



... 1 r 

P + 2HP - AitGqhP - — dr a 2 (r)TZ(X, r) 

a J to 



1 

ag 



-v, -n 



(47) 



For details on the possibility and advantage of this choice, 
see Adler & Buchert 1999, Appendix A. 



The left-hand-side of this equation reproduces the 'dust 
model'. If the correction to 'dust' is expressed in terms of 
the fields g, u, w, then the equation can be written entirely 
in terms of P and TZ by virtue of Eqs. (|39l46l51j) . Thus, 
the difficulty in deriving a closed, local differential equa- 
tion for the deviation field P lies in the term TZ, which is 
a nonlocal (in space and time) and a non-linear functional 
of P given by Eq. ifTOb l . 

As a further illustration of the general case we write 
down the corresponding equation for the density. For this 
purpose we introduce the contrast density field A := 
(Qh — q)/q, —1 < A < oo, which is more adapted to the 
non-linear situation than the conventional density con- 
trast 8 = (g — gn)/gH (for the 'dust' case see: Buchert 
1989, Eq. (7)ff, 1992 Eq. (31) ff, 1996 Eq. (25) ff (with a dif- 
ferent sign convention for A)). The exact evolution equa- 
tion for this variable can be obtained from Eqs. (|13b .b'l: 



A + 2HA - 4irGg H A = 



(1 + A) 



2IP C 



1. 



ga 



v g -n 



, (48) 



with the second scalar invariant of the peculiar-velocity 
gradient, 

//P0C ' = 2~0? ^ Ui ^ 2 ~ Ut '3 U i^ = ^Vq-(ttVq-M-M-VqM) . 

For the purpose of evaluating the Eulerian spatial 
derivatives in the right-hand-side of the above general 
equations, we employ the following transformation: 



0_ 



d 



where J 1 is the inverse to the following matrix: 

, dPi(X,t) 



Jij (X, t) 



dXi 



(49) 



(50) 



In the rest of this section we assume always that the 
map X — > q is invertible, i.e., no shell-crossing, so that 
det Jij > 0. Mass conservation, Eq. (|13b L relates the den- 
sity field to P via the matrix Jij 9 : 



1 /"* 

w(X,t)=4nGag s P(X,t) + ^ dr a 2 (t)K(X , t) , g(X,t) 

a J t 



Q[XM) det J l3 (X,t) 
g H (t)/ det Jij (X,t), 



(51) 



where we take initial quasi-homogeneity for simplicity 
(g(X,t ) w i?fr(to), Jij(X,t ) w <Sy). 

6.2. Exact evolution equations II: plane symmetry 

The problem is simplified if one considers particular geo- 
metric settings with high symmetry. Although this leads 
to not fully realistic models, it can provide hints of the 
mathematical and physical properties of the models in the 
non-perturbative regime. Plane-symmetric models offer 

9 Of course, g solves Eq. igHJ, if P solves Eq. B71 . 
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an excellent test case for numerical simulations of multi- 
stream systems concerning a variety of properties includ- 
ing spatial scaling (e.g. Doroshkevich et al. 1980; Gouda 
& Nakamura 1988, 1989; Yano & Gouda 1998; Fanelli 
& Aurell 2002; Aurell et al. 2003; Yano et al. 2003; the 
assumption of spherical symmetry plays a complemen- 
tary role, e.g. Padmanabhan 1996.) In such configura- 
tions, the fields vary spatially along a single direction, and 
so render the Lagrangian approach particularly suitable 
for an extension into the non-perturbative regime. Plane- 
symmetric models simplify the problem enormously, since 
the residual term TL vanishes, and the 'dust' evolution 
remains in the Lagrangian linear regime at all times (im- 
plying in particular that the parallelism relation (|24() - 
stating a spatially constant factor of proportionality - 
holds asymptotically at large times). 

In a plane-symmetric configuration, the fields vary 
only along a single direction, say i = 1. Eqs. (|46I47I49I51|I 
then are simplified to (we hereafter drop the index i = 1): 



g(X,t) = g H (t) 



dP(X, t) 



d_ 

dq 

w(X,t) 



1 



dP(X,t) 



dX 

n -1 



8X 

ATrGag H P(X,t) , 

<9n 



d 

dX 



(52) 



P + 2HP-A-kGq h P = —F y — — — 

ag a gH OX 



Therefore, depending on the specific model, we get differ- 
ent equations: 



'Dust model': 



P + 2HP - 4nGg H P = 



(53a) 



'EJN model': inserting Eqs. 116f) particularized to the 
plane-symmetric configuration, one obtains: 



P + 2HP-AnGg H P 



1 d 



a 2 g H dX P \l + dP/dX J ' 
(53b) 

'SSE model': using Eqs. (IS). 



Qh \ 



2 9 {dP/dXf 



P + 2HP-,nGg H P = -BL 2 - p/ , 



4irGg H 



2(l + dP/dX) 2 1 + dP/dX 



(53c) 



Eqs. (|53|l are closed partial differential equations for the 
function P(X, t). Notice that the equation for the 'dust 
model' is linear, so that Lagrangian linearization already 
yields the exact solution. The equation for P(X, t) in the 
'EJN model' becomes also a linear equation for the special 
dependence p cx g^ 1 (Chaplygin gas, App.|0 c/. Adler & 
Buchert 1999). This case is formally of interest because it 
provides a subcase in which the Lagrangian linear solution 
exactly solves the plane-symmetric problem. 



Eqs. should allow one to explore the role of the 
corrections to 'dust' beyond the assumptions employed in 
^3.11 in particular, the back-reaction on the 'dust' tra- 
jectories by the correction terms. It should be noted that 
Eqs. I|53|) are of second order in time, in contrast to the 
generalized 'adhesion model', Eq. I|26fa ,') — the relevance 
of this fact concerning the formation of singularities is dis- 
cussed in M6.4I In this connection Eqs. Q53b .c) have fea- 
tures of a quasilinear hyperbolic ("wave") equation, while 
Eq. (j2Ht) is a quasilinear parabolic ( "diffusion" ) equation. 
Gotz (1988) considered the 'EJN model' with the "isother- 
mal" relationship p cx g in the absence of background cos- 
mological expansion (a(t) = 1), and he was able to show 
that it can be mapped exactly to the Sine-Gordon equa- 
tion, which is known to admit soliton solutions. 

Fanelli & Aurell (2002) studied numerically the one- 
dimensional Vlasov-Poisson system (see |J2J j in particular 
the time-dependence of a collapsing isolated perturbation. 
Fanelli & Aurell tried to interpret the numerical findings 
in the framework of the simplified 'EJN model' l|26(l with a 
polytropic relationship for the pressure, p(g) cx g 1 . A sim- 
ple scaling argument led them to the value 7 = 7/9 (the 
original 'adhesion model' corresponding to 7 = 2). This 
conclusion is, however, questionable because boundary- 
layer theory cannot be applied to Eqs. (|2"BT when 7 < 1. 
A more general analysis (Domfnguez, unpublished) shows 
that a shock can form in this case only when the initial 
velocity is smaller than the maximum velocity, and this 
maximum velocity vanishes in the limit fi — > 0. That is, 
when 7 < 1, the adhesive term is too weak to prevent sin- 
gularities for arbitrary initial velocities. Thus, it might be 
the case that Eqs. (|2fi[) cannot describe even qualitatively 
the numerical experiment by Fanelli & Aurell. 

6.3. A new non-perturbative approximation 

As remarked, one of the difficulties with the exact Eq. (|4Tp 
is the nonlocal nature introduced by TZ. A possible ap- 
proximation consists in setting 1Z — 0, so that Eq. I|47|l 
simplifies to 



P + 2HP - 4nGg H P = 



og 



-Vq-n 



(54) 



As long as F and II are modelled as local functions of 
g, u and w (as in the models we consider in this work), 
this is a partial differential equation for P by virtue of 
Eqs. (|39I46I51II . albeit a difficult one. 

One can view the approximation 7Z. = in the same 
spirit as the parallelism assumption (|24|l : one gets rid of 
the problem of nonlocality by approximating the exact 
peculiar-gravitational field by a relationship which holds 
in the (both Eulerian and Lagrangian) linear regime of 
the 'dust model' as well as exactly in the plane-symmetric 
configuration: 

w = AirGag H P, (55) 

to be compared with w = 4nGagH{b/b)P correspond- 
ing to Eq. I|24f) . The relation 155fl replaces the "slaving" 
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assumption of parallelism of w to the peculiar-velocity. 
It is more general and simpler, because it does not in- 
volve specification of the function b(t). Unlike the con- 
dition of parallelism, this approximation also retains fea- 
tures of the exact field equations for w: e.g., at any given 
instant of time, the approximated field w depends only 
on the displacement (i.e., the distribution of mass) at this 
time, and not on its time-derivative (the peculiar-velocity 
field). This renders the differential equation obeyed by P 
second-order in time (see flfi.41 in this respect). A proper 
understanding of the extrapolation of the linear relation- 
ship between w and P into the non-perturbative regime 
requires that the residual term 1Z. be analyzed, a task be- 
yond the goals of this work. 

To get a glimpse of the mathematical difficulties posed 
by Eq. (|51Jl . we particularize it to the 'EJN model' of F 
and II, Eqs. Q16[l. The right-hand-side can be written, in 
analogy with Eq. (|17b 'l. as: 



1 

ag 



F--v q n 

a 



^A qW = gH { ig) A q P. 



2 

a z g 



(56) 



(In the expressions for F and II, the fields g and w are in- 
terchangeable via Eqs. (fHfc .dL but the resulting different 
equations are no longer exactly equivalent when the ap- 
proximation 1Z- — is used.) This term can be simplified 
with the "virial equation of state" p = ng 2 , so that the 
prefactor of the Laplacian is density-independent. This 
expression for p(g) is certainly a good phenomenological 
model consistent with the results of numerical simulations 
(Dommguez 2003, Dommguez & Melott 2004). The re- 
sults presented in H3. II also suggest that the precise value 
of the polytropic exponent does not alter the qualitative 
behavior in high-density regions. With this simplification, 
Eq. Q54JI becomes: 

P + 2HP - A-kGqhP = ^f^AqP. (57) 

This equation, being non linear in Eulcrian and in 
Lagrangian space, belongs to the class of hybrid 
Lagrangian/Eulerian models mentioned in Hfi.ll As such, 
the approximation may be viewed — in contrast to 
the Lagrangian perturbation scheme — in closer cor- 
respondence to the former 'adhesion approximation'. 
By construction, this approximation reduces to the full 
Lagrangian linear part of the 'dust model', Eq. I|42l) . 
when the correction to 'dust' is dropped. Upon linearizing 
the Eulcrian Laplacian, we recover the Lagrangian lin- 
ear approximation, Eq. (|45H . (In other words, our non- 
perturbative approximation extrapolates the Lagrangian 
linear model by the replacement Ax — * Ax+p.) 
Linearization of the convective non-linearities yields the 
Eulcrian linear approximation, Eq. (|34|1 (for this, note that 
linearization of Eq. (|51|l gives S = — Vo ■ P). 

In conclusion, we propose the approximation i|55|) as 
a first extrapolation into the non-perturbative regime of 
the Eulerian and Lagrangian linear approximations with- 
out resort to the constraining assumption of parallelism. It 



should be viewed as a simplification of the exact Eqs. I|13fl 
in the sense that the approximation yields a local equa- 
tion, which should facilitate the theoretical analysis. An 
attempt to solve the approximate equation could be suc- 
cessful, but lies beyond the scope of the present work. 

6.4. Regularization of singularities 

A point one would like to be able to prove rigorously be- 
yond the argument offered in HH.ll is whether the 'EJN and 
SSE models' do indeed prevent the formation of singu- 
larities given initial conditions of cosmological relevance. 
The absence of singularities when initially smooth data 
are propagated by the Vlasov-Poisson system has been 
proved mathematically (Lions & Perthame 1991; Schaeffer 
1991; Pfaffelmoser 1992; Rein & Rendall 1994). The origi- 
nal 'adhesion model' reduces to the 3D Burgers' equation 
(Eq. H26fl with constant /i), which has the unusual prop- 
erty of being solvable analytically. From there, one can 
demonstrate that no singularity arises for any sufficiently 
smooth initial conditions (see, e.g., Frisch and Bee 2001). 
For other models this is a very difficult question (e.g., it is 
still an open problem for the 3D incompressible Navier- 
Stokes equations (Frisch 1995, Sec. 9.3)), and here we can 
only provide some general remarks. 

Starting from Eqs. (|13l I21JI . one can derive the follow- 
ing evolution equations for the density g and the peculiar- 
expansion rate 9: 



g = —3Hg 



(58a) 



-2H9 - i<9 2 - o-ijo-ji + \y 2 - \itGa(g - g H ) 
+ -V„- ( -F — V„ • IT 



q \ v q 

a \g ag 



(58b) 



The continuity equation (|58h ,^ can be integrated in 
Lagrangian coordinates: 



g{X,t) = g(X,t ) 



exp 



to a ( T ) 



(59) 



Starting from bounded initial conditions, a singularity 
in the density field can occur only when 9 — ► — oo at 
some finite time 10 . From Eq. H58b i we see that only three 
terms can actually prevent 9 from becoming too negative: 
the background expansion, the vorticity, and perhaps the 
model-dependent correction to 'dust'. In particular, for 
the 'dust model' with vanishing initial vorticity, only the 
term linear in 9 in Eq. (|58b % ) can oppose collapse, so that 
a density singularity is unavoidable provided there is an 
initial overdensity. 

For the 'EJN model', the correction to 'dust' reads 
— Vq ■ (g _1 Vqp), where p is modelled as a density- 
dependent pressure, so that Eqs. I|17|l have the structure of 

10 In this work we are not concerned with more general (non- 
Lagrangian) singularities in the velocity-gradient field, such as 
those that may occur in incompressible flows. 
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the equations of inviscid fluid dynamics. In the absence of 
self-gravity, it is known that the solution becomes non- 
differentiable for a generic class of smooth initial condi- 
tions, so that the Lagrangian-to-Eulerian map q(X,t) 
ceases to be defined. The map can remain uni- valued if 
we allow for shocks (moving discontinuities in the fields): 
this is called a weak solution, and they have to be under- 
stood as solutions of the differential equations (|17|l with 
an additional viscous term fiA q u in the limit fi — > + , in 
the same way as we discussed for the adhesion model 11 . 
The presence of self-gravity is not expected to alter this 
conclusion; the question is rather whether the initial con- 
ditions of cosmological interest belong to the class of ini- 
tial conditions inducing shock formation. The reason for 
the emergence of shocks can be ultimately traced back to 
the fact that density and peculiar-velocity are indepen- 
dent variables: when convection by the peculiar-velocity 
field occurs at such a high rate that the density gradi- 
ents are not enough to build up a counteracting pressure, 
a singularity arises. However, in the cosmological context 
the peculiar-velocity field becomes "slaved" to the den- 
sity field in the linear regime (see ^4.1|l . The assumption 
of parallelism (12411 extrapolates this "slaving" relationship 
into the non-linear regime: this allows one to express the 
density-dependent correction to 'dust' in terms of the ve- 
locity gradient, and it seems sufficient to avoid singular- 
ities (see M3.1fl . But we have seen that a departure from 
the condition of parallelism is generated during the evolu- 
tion: an open question is therefore whether the assumption 
of parallelism is qualitatively good or instead the dynami- 
cally generated departure from it can lead to the formation 
of singularities. 

The mathematical properties of the equations underly- 
ing the 'SSE model' have been barely studied in the litera- 
ture. The correction to 'dust' is a complicated function of 
the fields. Unlike the 'EJN model', the corrections depend 
explicitly on the gradients of the peculiar-velocity field 
too, and the analogy with shock formation in inviscid flu- 
ids is not valid. Another important difference to the 'dust' 
and the 'EJN' models is that vorticity is generated by the 
'SSE' corrections to 'dust' (see H4.2|) . which, as pointed 
out, tends to oppose collapse. Vreman et al. (1996, 1997) 
have studied numerically the hydrodynamical equations 
of an ideal gas with an additional stress tensor given by 
Eq. I|19b ') when the initial condition is a perturbation to 
a smooth profile. An instability was discovered for which 
the 'SSE'-term is made responsible. Vreman et al. (1996) 
also undertook a theoretical analysis of the following sim- 
plified model (one-dimensional forced Burgers equation 
with a simplified, 'SSE'-like correction; ry > 0): 



du du d 2 u d ( du\ 



(60) 



11 The weak solutions can be interpreted equivalently as so- 
lutions of the integral form of the balance equations for mass 
and momentum. 



It is found that a sinusoidal profile (which fixes the form 
of the forcing f(x) by demanding stationarity) is linearly 
unstable to perturbations if 2r] > fJL—1. Indeed, as the dis- 
cussion of Eq. l|6"2|l below exemplifies, it is not true that 
the 'SSE'-corrections will behave like viscous terms un- 
der all conditions. As with the 'EJN model', the ques- 
tion is whether the 'SSE model' does indeed prevent the 
formation of singularities when the initial conditions are 
restricted to be of cosmological relevance. 

An issue closely related to this problem is the behavior 
of the correction to 'dust' concerning the local "energy 
budget" . The density of mean peculiar-kinetic energy is 
defined as K :— (l/2)g\u\ 2 . From Eqs. ifHfl) one finds: 

9if 1„ , „, 

(61) 



dt 



-V 



q ■ (uK) + 5HK = u-(gw + F) 



1 



1. 



n : (y qU ) - -v q • (n • u) , 



a a 
where : indicates contraction of two indices. The first term 
on the right-hand-side represents the work done by the 
total (mean field and corrections thereof) gravitational 
force, the second term is the work done by the velocity dis- 
persion, and the third term is the contribution due to the 
flux of velocity dispersion. In the 'EJN model', Eqs. Ijl6|) . 
only velocity dispersion does work: II : (Vqii.) = p(g)0. 
Given that p > 0, this term tries to reduce the kinetic en- 
ergy in (relative to the Hubble expansion) locally collaps- 
ing regions (8 < 0) and to increase it in locally expanding 
ones (8 > 0). 

In the 'SSE model', the sign of the work done by the 
deviations from mean field gravity, u ■ F = BL 2 (uV ' q g) : 
(W q w), depends on the orientation of the tidal tensor 
and, like the work done by mean field gravity, is in gen- 
eral unrelated to whether the fluid is locally compress- 
ing or expanding. On the other hand, the work done by 
velocity dispersion is most easily expressed in the basis 
that diagonalizes the symmetric part of (l/a)diUj (this is 
(Tij + (l/3)85ij, according to the decomposition J2U): 



In 

a 



-BL g(d k u i )(d k u j )(d i Uj 



(VqM 

= 2B(aL) 2 gY J S^ 



'(0 



(62) 



where are the eigenvalues of the tensor cry + (1/3)9 Sij , 
satisfying S(j) = 8. The work is the sum of contribu- 
tions along the three eigendirections: the contribution of 
the direction of local compression (s^ < 0) tends to re- 
duce the kinetic energy, that of the direction of local ex- 
pansion (s(j) > 0) tends to increase it. The net sign is the 
outcome of this competition. In this sense, the 'SSE' ex- 
pression for the velocity dispersion can be viewed as a "di- 
rectional" generalization of the case of an isotropic pres- 
sure like in the 'EJN model'. It must be noticed, however, 
that this competing effect between local expansion and 
compression is effective only in the fully three-dimensional 
case: it can be easily checked that in the restricted cases 
of one or two-dimensional motion, expression 

(E3 has the 

same sign as 9. 
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7. Summary and prospects 

Starting from the equations of motion for point parti- 
cles, we have implemented a smoothing procedure in phase 
space leading to equations for the smoothed density and 
peculiar-velocity fields in Eulerian space. These equations 
are the first ones in an infinite hierarchy of equations, and 
approximations to close the hierarchy are required. A clo- 
sure approximation corresponds physically to an ansatz 
on how the evolution of the filtered fields is coupled to the 
degrees of freedom below the filtering scales. In this uni- 
fied context, we have addressed three different closures, 
leading to three different models proposed earlier in the 
literature: the 'dust model', the 'EJN model', and the 
'SSE model'. The 'dust model' describes the evolution of 
fluid volume elements under the sole influence of the self- 
consistently generated gravitational field. The corrections 
to 'dust' take into account the dynamical influence of the 
small-scale degrees of freedom. The 'EJN model' approxi- 
mates this influence as a density-dependent pressure. The 
'SSE model' expresses this influence as a functional of den- 
sity and velocity in a systematic expansion. 

We have provided a clear explanation of how the cor- 
rections to the 'dust model' by both the 'EJN closure' 
and the 'SSE closure' lead to models that are qualita- 
tively equivalent to the 'adhesion model'. The key hy- 
pothesis is the parallelism relation (|24|l between peculiar- 
velocity and peculiar-acceleration. This already holds in 
Zel'dovich's approximation to the exact 'dust' evolution — 
the assumption underlying the 'adhesion approximation' 
consists of extrapolating it to also evaluate the corrections 
to 'dust'. 

We have shown that within both the Eulerian and 
Lagrangian perturbation frameworks, the corrections to 
'dust' violate this assumption in general. Going beyond 
the phenomenology of the 'adhesion model' therefore re- 
quires relaxing the parallelism assumption. When this is 
done, differences between the 'EJN model' and the 'SSE 
model' appear: deviations from the parallelism relation- 
ship arise at linear (both Eulerian and Lagrangian) per- 
turbative order in the 'EJN model', but not in the 'SSE 
model'; vorticity is generated by the 'SSE'-corrections to 
'dust', but not by the 'EJN'-corrections. 

The non-perturbative theoretical investigation raises 
considerable difficulties. We suggested the plane- 
symmetric case as the "exact body" of the non- 
perturbative equations that may serve as a starting point 
for a deeper understanding of the presented models. This 
case can be reduced to the study of a single partial dif- 
ferential equation for the displacement field, offering the 
possibility to study the 'EJN' and 'SSE' models beyond 
any approximation. This simplification misses some fea- 
tures that are expected to be of relevance for the de- 
tailed inner structure of collapsing high-density regions 
(e.g., generation of vorticity or the fact that velocity and 
gravity have different directions). We suggested a new 
non-perturbative approximation, which yields a hybrid 
Eulcrian/Lagrangian partial differential equation for the 



displacement field that contains as limits (i) the Eulerian 
and Lagrangian linear approximations, and (ii) the stan- 
dard 'adhesion approximation' if one imposes the assump- 
tion of parallelism of peculiar-velocity and -acceleration. 
When particularized to the 'EJN model', the 'virial equa- 
tion of state' p oc q 2 provides substantial simplifications 
and, at the same time, is a reasonable choice. We argued 
how this approximation can be viewed as the next step 
beyond the parallelism assumption required by the 'adhe- 
sion model'. Finally, we provided a qualitative discussion 
of the exact 'EJN' and 'SSE' models concerning the for- 
mation of singularities. 

An issue we have not addressed in this work is the 
choice of the smoothing scale L. This must be guided by 
the expected applicability of the corresponding closure ap- 
proximations, in turn likely related to the initial condi- 
tions. The investigation of this point seems to require first 
the ability to extract predictions from the new models in 
order to validate them with respect to observations and 
viable N-body simulations. 

Further exploration and numerical realization of the 
models presented will hopefully lead us not only to ap- 
proaches capable of following large-scale structure into 
the non-linear regime, but will also offer an interesting 
way of understanding the role of previrialization (Peebles 
k Groth 1976, Davis & Peebles 1977, Peebles 1990, Lokas 
et al. 1996), and the onset of virialization in connection 
with galaxy and galaxy cluster formation. Furthermore 
such models can be employed to improve generic inhomo- 
geneous collapse models (Buchert et al. 2000; Kerscher et 
al. 2001), and to enhance the power of reconstructions 
of initial data from present-day density and peculiar- 
velocity fields (e.g., Croft & Gaztahaga 1997; Brenier et al. 
2003; Mohayaee et al. 2003 for recent efforts). In the latter 
context, the proposed non-perturbative approximation, 
based on the relation w oc P, Eq. (55} ) ma y provide an 
alternative to the assumption of parallelism of peculiar- 
acceleration and peculiar-velocity, w oc u, Eq. 1)24(1. 
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Appendix A: On the transformation to comoving 
coordinates 

To derive Eqs. (|13f) we first introduced comoving coordi- 
nates, Eqs. i|10fl . and then the coarse-graining procedure 
was applied. It is legitimate to ask whether the order can 
be exchanged, i.e., what are the transformations that lead 
from Eqs. JHJ to Eqs. (jl3(l ? They can be obtained easily by 
inserting the transformations (|10l) into the definitions l|12l 
I14|l and using the definitions {7| |§J). One gets (here we 
write explicitly the dependence on the smoothing length) : 



1 

-x, 

a 



L = -C , 
a 



g(q,t;L) = p(x,t;C) , 



u(q, i; L) = v(x, i; C) — HR(x, t; C) 



(A.la) 



(A. lb) 



(A.lc) 



w{q,t; L) = g(x,t; L) H x . (A. Id) 



F(q, t; L) = ?{x, t; C) + 4nGg " — -{R - x)p{x,t; £) , 



(A.le) 



Ihj (q, t; L) = Vij (x, t; C) - H[V Z] (x, t; C) + V 3l (as, t; £)] 



-H 2 l i:j {x,t;C) 



L., i \ .1 . I . I . (A. If) 

The following new quantities had to be introduced: 

„n m v—v fx — x^\ /„•, 
pR(x,t;jC):=-^J2 W \ r ) * ' ( A ' 2a ) 



C 



2^. (*,*;£) :=^J2 W 



-pviRj(x,t;C), (A.2b) 
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2y (*,*;£) :-^E^ 



X X^ ^ \ ( a \ ( a \ 



-pRiRj{x,t;C). (A.2c) 



i? is the center-of-mass position of the subsystem defined 
by the smoothing window, and X is the inertia tensor of 
this same subsystem with respect to R. Concerning 27, we 
note that the antisymmetrized combination tijkDkj is the 
angular momentum of this subsystem with respect to R; 
the transformation (|A.lf ) involves only the symmetrized 
tensor. Therefore, to obtain a self-contained set of trans- 
formations, new dynamical equations for R, X, T> would 
be required, including also further approximations to close 
the corresponding hierarchy. The simplest one is provided 
by the 'dust model': the small-scale degrees of freedom 
are neglected altogether as irrelevant for the dynamics, so 
that R — x, X ~ 0, T> — 0, and the transformation (|A.1(I 
reduces to the standard one (e.g., Peebles 1980). 

Compared to this 'standard transformation', transfor- 
mation (|A. 1|) differs in two points: (i) the additional trans- 
formation rules IjA.lb .f) involving the fields F and n, 
and (ii) the transformation between u and v, because 
in general R(x,t) ^ x. When the particle distribution 
looks homogeneous at all cosmological scales, the devia- 
tion R — x is irrelevant. In the highly non-linear regime 
of large inhomogeneities, one expects \R — x\ ~ aL, but 
then ^S> HaL nevertheless, and the cosmological 

expansion is also irrelevant (Dommguez & Gaite 2001). 
The difference may be important in the most interesting 
case when the structure is beginning to enter the non- 
linear regime, \v\, \ u\ ~ HaL. 

Although the change to comoving coordinates does not 
affect the physical content of the resulting equivalent equa- 
tions, © or (|T3)l . this equivalence can be broken after ap- 
proximations are made. It seems more straightforward to 
"subtract" the known background cosmological expansion 
before performing any smoothing/averaging, as we did to 
derive Eqs. lfT3)l. so that the approximations are concerned 
only with the "unknown" part of the evolution of the 
inhomogeneities. However, the above remarks show that 
this procedure "hides" additional physics and requires a 
careful definition of "background" at the coarse-graining 
scale, which is only implicit by coarse-graining in the pre- 
defined comoving space. 



Appendix B: Derivation of the 
'Small-Size-Expansion (SSE)' 

In this Appendix we recall the derivation of the 'SSE' 
closure (|19|) (Dommguez 2000). In analogy with defini- 
tions H12J) . one can formally define microscopic fields by 
integrating the Klimontovich density ® over the velocity 



variable 12 : 
Qmic(q,t) : 



mJduf K {aq,u-a q ,t) 



m 

773 



(B.la) 



m J du u fx(aq,u — aq,t) 
^E« (o) %-9 (a) ), ( B - lb ) 



W r 



Jin if y *' in a 



| *)/*(«* «-«**) 
-E™ (a)< 5(<?-g (a) ) . ( B - lc ) 



already expressed in comoving coordinates. In Eq. l|B.lc|) . 
w(q,t) := g(aq,t) + (4nGgff(t) — A)aq/3, with g given 
by Eqs. I@p). Note that w m i C (q,t) = w(q,t), so that by 
virtue of Eqs. fllfr) (when expressed in comoving coordi- 
nates), we can write: 



—4irGa(g n 



Qh) 



Vq x w mlc = 



(B.2) 

In terms of the microscopic fields, definitions (|12II14|I can 

be rewritten as follows: 



gu(q,t) 



dq' 



W 



L 



Qvnic 

(q',t), (B.3a) 



(q',t), (B.3b) 



(q',t) - gw(q,t) , 
(B.3c) 



U(q,t) = J JJ W ^ ) QmicUmicUmic(q' ,t) 

— guu(q, t) , (B.3d) 

where w is the mean field strength of the peculiar-gravity, 
defined by Eqs. (|13b .d). The introduction of the mic-fields 
is a formal step to express the definitions in terms of vol- 
ume integrals over fields. Notice in particular that u is 
defined by a volume average of the peculiar-momentum 
density (i.e., by a mass-average of the velocity) and thus 
it has the physical meaning of center-of-mass velocity of 
the coarsening cell. These fields obey Eqs. (|B.3b . b), which 
can be formally inverted to yield the following expansion 
in L (this is most easily derived in Fourier space; this for- 
mal inversion is defined provided the Fourier transform of 
the window W(-) has no zeros): 



gmic 



l- l -B{LV q Y + 0{LV q ) 4 



Q 



(B.4a) 



12 As with Eqs. iTJJ, the formal definitions of the fields 
Umic(q), Wmic{q) can be regularized by smearing Dirac's delta 
function. 
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Special values of "7 



with 



1-\b{L V q ) 2 - BL 2 ^ ■ V q + 0(L V q ) 4 



U 



JdzzW{z)=0, B := ^ y dz z 2 W(z) 



(B.4b) 



(B.5) 



These expansions show how the microscopic fields can be 
recovered by taking into account finer and finer details 
in the spatial distribution of the coarse-grained fields. To 
derive Eq. (jl9b ). one inserts Eqs. (|B.4|I into the defini- 
tion (|B.3f O encompassing the coupling to the small-scale 
kinetic degrees of freedom: 



U(q,t) 



dq^_ 

L 3 



W 



{guu(q' , t) — guu(q, t) 



- -BL 2 [uW 2 q ,(gu)(q',t) + {V 2 q ,u)gu{q' ,t) 

+ 2{V q ,g) ■ {V q ,u)u{q\ t)]} + 0(L 4 ) 
= BL 2 g{d k u)(d k u)(q 7 t) + 0(i 4 ) , 

where the fields in the integrand have been expanded 
about q 1 = q consistently to a given order in L. The ex- 
pansion l|19fa) for the correction to mean field is derived 
in the same manner (for this, notice that w m i C is a linear 
functional of g m i C given by Eqs. (|B.2|l h 



Appendix C: Dynamical equations of state and 
corresponding model features 

An ingredient of the 'EJN' model is an "equation of 
state" 13 p(g) for the trace of the velocity dispersion tensor 
as function of the mass density, Eq. I|16ll . In this Appendix 
we summarize the evidence in favor of this assumption and 
the physical meaning of some particular functional depen- 
dences. 



Numerical simulations 

The density-dependence of the trace of the velocity dis- 
persion i|14b) was addressed systematically with N-body 
simulations of CDM scenarios by Dommguez (2003) and 
Dommguez & Melott (2004) and can be summarized in 
a (sometimes piecewise) polytropic dependence p oc g 7 , 
which should be understood as an average, because the 
numerical data scatter around this dependence. The in- 
dex 7 is always found to lie in the range 1 < 7 < 2. The 
more power there is in the small scales of the initial mass 
distribution relative to the large scales, the smaller is 7. 



13 This notion as well as the names attached to the different 
polytropic models (with the polytropic index n defined through 
7 =: (n + l)/n) is borrowed from thermodynamics. We empha- 
size that here we refer to dynamical stresses and do not imply 
thermal equilibrium. 



There are values of the index 7 which have a particular 
significance (see also Buchert et al. 1999): 

• 7 = 2: 'Virialized state' 

This model is frequently encountered in our paper. From 
a mathematical point of view, the index 7 = 2 entails a 
simplification because the length X(g), Eq. 118JI . becomes 
p-independent. Thus, the standard 'adhesion model' is 
recovered in this case because the GM coefficient /1, 
Eq. J26bj), is also ^-independent; see also the derivation 
of Eq. (|57|). 

Physically, one can argue on the basis of the virial theo- 
rem (e.g. Chandrasekhar & Lee, 1968): a well-known con- 
sequence of it is that, if a system dominated by gravity 
is isolated and in a stationary state, then the total kinetic 
energy IA (in the frame of the system's center of mass) and 
the total (gravitational) potential energy W are related by 
the equation HA + W = 0. If our system is a coarsening 
cell of size ~ L defined by the smoothing window W(-), 
we have IA oc p and W ~ GL 5 g 2 assuming that W can be 
estimated reliably within a mean field approximation (i.e., 
not too much small-scale structure). To the extent that a 
coarsening cell can be considered approximately isolated 
and stationary, the virial theorem implies p oc g 2 . The 
use of this relationship in the evolution equation (|17b ) al- 
ways supposes an approximation, because the term — X7p 
represents the flux of momentum by the exchange of parti- 
cles between neighboring coarsening cells, in contradiction 
with the assumption of isolation required to derive the de- 
pendence p oc g 2 . 

• 7=1: 'Isothermal state' 

In this case, the plane-symmetric problem of the 'EJN 
model', in the absence of cosmological expansion, can 
be mapped exactly to the 'Sine-Gordon equation' (Gotz, 
1988). 

Physically, this dependence implies that the kinetic en- 
ergy per particle, which is proportional to pj g, is density- 
independent. In the context of a fluid in thermal equilib- 
rium, the reason is that, concerning any single particle, 
the rest of the system behaves like a thermostat. 

• 7 = 5/3: 'Adiabatic state' 

Starting from the definition l|14b h one can compute the 
following evolution equation for the tensor Iljj : 

an 



--U ik {d k Uj) - -Il jk (d k Ui) (C.l) 
a a 

— d k {u k Uij) dkCijk + Vij. 

a a 

In this equation, Cij k is the third reduced velocity moment 
(representing the heat flux in dilute gases) and "Py involves 
the departures from mean field. (The precise expressions 
for these two tensors can be found in (Dommguez 2000)). 
One may take the closure assumption that these two ten- 
sors are negligible, meaning physically that corrections to 
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mean field are unimportant, and that velocity dispersion 
is relatively small (and with it the "heat flux" =>■ "adia- 
batic evolution" ) . If one further makes the isotropy ansatz 
ILij = pdij, this equation in combination with the conti- 
nuity equation l|13b) then yields p = ng 5 ^ 3 and the strong 
kinematical constraint of vanishing shear, diUj+djUi oc 5ij 
(Buchert & Dommguez 1998). Note that this implies that 
the only consistent exact solution is spatially homogeneous 
(Triimper 1967; Ehlers & Rienstra 1969; the homogeneous 
model also satisfies 'adiabaticity'.) 

The coefficient of proportionality K is a local function of 
the Lagrangian coordinates given by the initial data; if we 
assume that, initially, this function is a spatial constant, 
then we are entitled to talk about a "state", since the 
relation between dynamical "pressure" and mass density 
is then global. 



(Dommguez, unpublished) demonstrates that "pressure" 
is too weak to avoid the occurrence of singularities when 
7 < 1. 



• 7 = 3: Plane symmetric 'adiabatic state' 
In a plane-symmetric configuration, Eq. I|C.1|) can be in- 
tegrated under the assumption that Cijk and Vij are neg- 
ligible without the need for the additional isotropy ansatz 
for Ily. Thus, the equation for the single component Tl xx 
becomes 

dU xx 1 dU xx 3^ du 

a . + ~u—^ h 5HU XX = — U xx — , 

at a ox a Ox 

from which one obtains Ti xx oc g 3 . 

Without regard to this derivation, in the context of the 
generalized 'adhesion models', the g-dependence of the 
GM coefficient derived with the 'SSE model', Eq. P5)l. 
can be obtained with the 'EJN model' and a dependence 
p oc g 3 , Eq. l|26bf . 



• 7 — 4/3: 'Cosmogenetic state' 

If one looks for a solution of Eqs. (|17fl in "comoving hy- 
drostatic equilibrium", i.e., u = but \7g ^ 0, with an 
"equation of state" p = ng 1 , then one can easily check 
that 7 = 4/3 is the only possible index for which a so- 
lution exists when n is time-independent (Chandrasekhar 
1967). 



• 7 = —1: 'Chaplygin state' 

In our context, this formal case is of interest because the 
particularization of the 'EJN model' to a plane-symmetric 
configuration gives a linear equation for the Lagrangian 
displacement field P exactly, see H6.2I 
More generally, the index 7 = — 1 introduces mathemat- 
ical simplifications in the integration of Eulcr's equation. 
Recently, the Chaplygin gas has received attention in cos- 
mological research as one possible model of dark energy 
(see e.g. the recent preprint by Gorini et al. 2004 and 
refs. therein). 



• 7 < 1: 

Application of boundary-layer theory to Eqs. (|26|) 
(Domfnguez 2000) shows that density singularities are reg- 
ularized when 7 > 1. On the other hand, a similar analysis 



